{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particle Swarm Neural Network Tutorial\n",
    "\n",
    "This is a tutorial to implement a feedforward neural network using particle swarm optimization (PSO) to train the network. The \"go-to\" approach for training a neural network is typically backpropagation, but there is merit to experimenting with metaheuristics to optimize neural network weights.\n",
    "\n",
    "I've previously written a tutorial on both <a href=\"https://www.github.com/stratzilla/neural-network-tutorial\">neural networks (using backpropagation)</a> and <a href=\"https://www.github.com/stratzilla/particle-swarm-optimization-tutorial\">particle swarm optimization</a> and want to combine both into one to show how PSO can be used to train a neural network in lieu of backpropgation.\n",
    "\n",
    "This tutorial does presuppose you've read both the linked tutorials above, so I won't be going into as much detail on how neural networks or PSO work alone. While I start my PSO-trained network from scratch, I don't restate things that I've said previously. Feel free to distribute, edit, make changes, improve, etc this tutorial. If you publish it elsewhere, please let me know and link back.\n",
    "\n",
    "## Neural Networks\n",
    "\n",
    "Without a complete retread of my neural network tutorial, you should know that for the purpose of this tutorial, our neural network is simply a series of weights between neurons. Unlike with backpropagation training, where the network stored neural error deltas and outputs alongside their weights, we only really care about the weights. However, the neural network doesn't optimize the synapse weights, rather we rely on the metaheuristic instead for this purpose.\n",
    "\n",
    "## Particle Swarms\n",
    "\n",
    "For our neural network, we'll be using PSO to train the network. How? To perform PSO on some optimization problem, we need to know the fitness function as well as the search space. For our neural network, the possibilities of synapse weight combinations is our search space and our error or loss function is our fitness function. We'll be using the mean squared error as our error function. Over time, PSO will find an optimized set of weights which we then initialize our network as.\n",
    "\n",
    "While we know what we're searching for within our search space, we don't necessarily know the bounds of the search space. We'll make an informed guess that our neural network weights will be fairly small and similar in magnitude so we can still perform PSO even without knowing the search space bounds.\n",
    "\n",
    "Like vanilla PSO, we know our search space (loosely) and our fitness function, so we can generate a swarm of candidate solutions to our optimization problem and then \"evolve\" them over time. Over time, solutions within our swarm will congregate around the global optima for our loss function.\n",
    "\n",
    "## Pseudocode\n",
    "\n",
    "Pseudocode for the algorithm is as below:\n",
    "\n",
    "```\n",
    "procedure pso-nn() is:\n",
    "  initialize network structure\n",
    "  initialize random swarm s of particles\n",
    "  for each epoch do:\n",
    "    for each particle p in swarm s do:\n",
    "      for each dimension d in search space do:\n",
    "        new_vel, new_pos = [], []\n",
    "        new_vel[d] := w * p.vel\n",
    "        new_vel[d] += c1 * r1 * (p.best_pos[d] - p.pos[d])\n",
    "        new_vel[d] += c2 * r2 * (s.best_pos[d] - p.pos[d])\n",
    "        new_pos[d] := p.pos[d] + new_vel[d]\n",
    "      endfor\n",
    "      p.vel = new_vel\n",
    "      p.pos = new_pos\n",
    "    endfor\n",
    "    for each particle p in swarm s do:\n",
    "      if fitness(p.pos) < p.fit do:\n",
    "        p.fit := fitness(p.pos)\n",
    "      endif\n",
    "    endfor\n",
    "    update swarm best\n",
    "    update particle bests\n",
    "  endfor\n",
    "  network weights := best from s\n",
    "endprocedure\n",
    "```\n",
    "\n",
    "Note this is identical to the pseudocode for the PSO algorithm as found in that tutorial, with the only addition being the best in swarm serves as weight initialization for the network.\n",
    "\n",
    "## Before you Begin\n",
    "\n",
    "This tutorial uses Python and two non-standard dependencies:\n",
    "\n",
    "- Python 3.6\n",
    "- GNU/Linux (or WSL, LXSS, etc)\n",
    "- X11 or equivalent\n",
    "- `matplotlib`\n",
    "- `pandas`\n",
    "- some data\n",
    "\n",
    "We'll be using `matplotlib` to graph the mean squared error and model accuracy over time but otherwise doesn't need this library. I use it for illustrative purposes but you can omit it if you don't wish to verify the system works. `pandas` is used for a very important step: normalizing the data. The algorithm will work without it but at a detriment.\n",
    "\n",
    "As for the data, I will be using the Wine data set (see `/data`) but the implementation will be data agnostic. Since the network's input and output layer structures are dependent on the data, it will be generalized to take any data set and adjust the model as necessary.\n",
    "\n",
    "## Loading the Data\n",
    "\n",
    "Before implementing the network and metaheuristic, we need to handle how data is loaded and stored.\n",
    "\n",
    "The Wine data set consists of 178 examples of 13-dimensional data (thirteen continuous features per row) along with a discrete classification (three classes). We'll split the data into training and testing sets with a 70%/30% ratio. This means the network is trained on 70% of the data and tested upon the remaining 30%.\n",
    "\n",
    "For our PSO algorithm to work, we need to know the search space. While we know the search space is the set of possible weight permutations, we don't yet know the boundaries of the search space. In fact, there's no declarative way to determine the space boundaries but we can make an informed guess that the weights will be relatively small and of similar magnitudes. We'll normalize our input data around 0 to shrink the total search space and make an assumption the weights will be within `[-7.00, 7.00]`.\n",
    "\n",
    "Here we'll load the data, normalize it, and then split it into our training and testing portions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "def load_data(filename):\n",
    "    \"\"\"Loads CSV for splitting into training and testing data.\n",
    "    \n",
    "    Parameters:\n",
    "        filename : the filename of the file to load.\n",
    "    \n",
    "    Returns:\n",
    "        Two lists, each corresponding to training and testing data.\n",
    "    \"\"\"\n",
    "    # load into pandas dataframe\n",
    "    df = pd.read_csv(filename, header=None, dtype=float)\n",
    "    # normalize the data\n",
    "    for features in range(len(df.columns)-1):\n",
    "        df[features] = (df[features] - df[features].mean())/df[features].std()\n",
    "    train = df.sample(frac=0.70).fillna(0.00) # get training portion\n",
    "    test = df.drop(train.index).fillna(0.00) # remainder testing portion\n",
    "    return train.values.tolist(), test.values.tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the above function, we load the CSV into a master dataframe, normalize the features, then take a random 70% sample as training data and use the remaining as testing data. We'll output the number of examples for each to verify as well as find the number of columns in the data (attributes) and the number of unique classifications in the data (classes):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training examples: 125\n",
      "Testing examples: 53\n",
      "\n",
      "Features: 13, Classes: 3\n"
     ]
    }
   ],
   "source": [
    "# get training/testing data\n",
    "TRAIN, TEST = load_data('./data/wine.csv')\n",
    "\n",
    "print(f'Training examples: {len(TRAIN)}')\n",
    "print(f'Testing examples: {len(TEST)}')\n",
    "\n",
    "# count number of attributes\n",
    "FEATURES = len(TRAIN[0][:-1])\n",
    "# count number of unique classifications\n",
    "CLASSES = len(list(set([c[-1] for c in (TRAIN+TEST)])))\n",
    "\n",
    "print(f'\\nFeatures: {FEATURES}, Classes: {CLASSES}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is congruent with what we know about the data set in terms of features and classes. We'll use these counts when initializing the network shortly as we want the input layer size to match the number of attributes and the output layer to match the number of unique classifications. The input and output layers are a consequence of the data and not user defined.\n",
    "\n",
    "## Initializing the Network\n",
    "\n",
    "Since we know the size of the input and output layers, we need only come up with a hidden layer size to begin network initialization. We'll use six hidden neurons in our lone hidden layer giving us a `13-6-3` network structure.\n",
    "\n",
    "Later on we'll need a way to encode a particle to use as weights for the network, so our network initializer will take this into consideration. We're still using bias for our network, so each of our particles will need to be `h * (n+1) + o * (h+1)` long, where `n` is the size of our input layer, `h` our hidden layer, and `o` our output layer. For our `13-6-3` network, this means our particles will be in `105` dimensions.\n",
    "\n",
    "Here is our network initializer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "HIDDEN_SIZE = 6\n",
    "\n",
    "def initialize_network(p):\n",
    "    \"\"\"Neural network initializer.\n",
    "    \n",
    "    Parameters:\n",
    "        p : the particle to encode into the network.\n",
    "    \n",
    "    Returns:\n",
    "        The n-h-o neural network.\n",
    "    \"\"\"\n",
    "    n, h, o = FEATURES, HIDDEN_SIZE, CLASSES\n",
    "    part = iter(p) # make iterator from p\n",
    "    neural_network = [] # initially an empty list\n",
    "    # there are (n * h) connections between input layer and hidden layer\n",
    "    neural_network.append([[next(part) for i in range(n+1)] for j in range(h)])\n",
    "    # there are (h * o) connections between hidden layer and output layer\n",
    "    neural_network.append([[next(part) for i in range(h+1)] for j in range(o)])\n",
    "    return neural_network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like how a network should be initialized with weights within a small interval such as `[-1.00, 1.00]` or `[0.00, 1.00]`, we want to also initialize our particles in `n`-space at random with each axis being within a similarly small interval. Right now we'll make a sample particle and output the network structure to verify it was initialized correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.68,  0.15,  0.60,  0.60,  0.18,  0.91, -0.50,  0.32,  0.79, -0.67,  0.04, -0.08,  0.19, -0.19]\n",
      "[ 0.49, -0.19, -0.63, -0.66,  0.31,  0.17, -0.56,  0.89,  0.18, -0.13, -0.93, -0.20,  0.25,  0.38]\n",
      "[ 0.13, -0.79, -0.35,  0.86, -0.04,  0.18, -0.51, -0.29,  0.58, -0.69, -0.97, -0.84,  0.81,  0.39]\n",
      "[-0.56,  0.51,  0.16, -0.07,  0.29,  0.79,  0.15, -0.77, -0.29,  0.73,  0.81, -0.56, -0.67,  0.37]\n",
      "[-0.95,  0.12, -0.79, -0.76, -0.02,  0.35, -0.90, -0.81,  0.04,  0.53, -0.02,  0.64, -0.53,  0.91]\n",
      "[ 0.36,  0.52, -0.98,  0.45, -0.65,  0.85,  0.19, -0.01,  0.27,  0.78, -0.63,  0.14,  0.95, -0.18]\n",
      "\n",
      "[-0.01,  0.30,  0.89, -0.01,  0.11, -0.27, -0.16]\n",
      "[ 0.32, -0.38,  0.65, -0.90, -0.67,  0.27,  0.45]\n",
      "[-0.60,  0.98,  0.35,  0.48, -0.29, -0.95,  0.60]\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "DIMENSIONS = (HIDDEN_SIZE * (FEATURES+1)) + \\\n",
    "    (CLASSES * (HIDDEN_SIZE+1))\n",
    "\n",
    "SAMPLE_PARTICLE = [round(random.uniform(-1.00, 1.00),2) \\\n",
    "    for _ in range(DIMENSIONS)]\n",
    "\n",
    "NETWORK = initialize_network(SAMPLE_PARTICLE)\n",
    "\n",
    "for layer in NETWORK:\n",
    "    for neuron in layer:\n",
    "        print(neuron)\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's been truncated for readability; when we begin making particles, we don't want to trim off excess decimal places.\n",
    "\n",
    "What you should take from this output is there are two layers of connections (between the input and hidden layers, and between the hidden and output layers). For the first set of weights, there is a set for each of the neurons in the hidden layer (`6`), each consisting of `14` weights (`13` from input layer + `1` bias). For the second set of weights, there is a set for each of the neurons in the output layer (`3`), each of `7` weights (`6` from hidden layer + `1` bias).\n",
    "\n",
    "This is congruent with our expected structure, however seeing the structure is not terribly important right now.\n",
    "\n",
    "## Fitness Function\n",
    "\n",
    "So far, we have a way to encode the candidate solution into a network but not a way to test how well the network performs. We've established our fitness function will be the mean squared error function, and to this end, we need to implement some of the busywork involved with any feedforward neural network.\n",
    "\n",
    "We'll need to feed data forward into the network (it is a feedforward network after all) to calculate the mean squared error, so much of the below is just an implementation of those requirements.\n",
    "\n",
    "### Summing and Activation Functions\n",
    "\n",
    "The summing function takes neural inputs and aggregates them to be sent to the next layer. In our network, for example, one of the hidden neurons takes eight inputs and needs to create one output from it. To this end, we use a summing function to concatenate them by using the weighted sum of the neural inputs. Each input will contribute a different amount to the output and we use the synapse weights to determine this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summing_function(weights, inputs):\n",
    "    \"\"\"Sums the synapse weights with inputs and bias.\n",
    "    \n",
    "    Parameters:\n",
    "        weights : synaptic weights.\n",
    "        inputs : a vector of inputs.\n",
    "    \n",
    "    Returns:\n",
    "        The aggregate of inputs times weights, plus bias.\n",
    "    \"\"\"\n",
    "    # bias is the final value in the weight vector\n",
    "    bias = weights[-1]\n",
    "    summ = 0.00 # to sum\n",
    "    for i in range(len(weights)-1):\n",
    "        # aggregate the weights with input values\n",
    "        summ += (weights[i] * float(inputs[i]))\n",
    "    return summ + bias"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have the weighted sum of the neural inputs, this sum is transferred to an activation function. This function determines the neuron activation or firing: for linear activations, this is a binary \"ON\" or \"OFF\", whereas with non-linear activations, it's a degree of \"ON\".\n",
    "\n",
    "We'll use the rectified linear unit activation function:\n",
    "\n",
    "$$\\text{ReLU}(z) = \\left\\{ \\begin{array}{ll} z & \\text{if }z \\geq 0 \\\\ 0 & \\text{if }z \\lt 0 \\end{array} \\right.$$\n",
    "\n",
    "We can implement this like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def activation_function(z):\n",
    "    \"\"\"ReLU activation function.\n",
    "    \n",
    "    Parameters:\n",
    "        z : summed output of neuron.\n",
    "    \n",
    "    Returns:\n",
    "        The neuron activation based on the summed output.\n",
    "    \"\"\"\n",
    "    return z if z >= 0 else 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I've chosen ReLU because some of the problems inherent with ReLU, such as vanishing gradient or dead neurons are a non-issue when we're not using the activation derivative. ReLU is a great activation function when these issues are non-existent, which makes it the perfect candidate for our activation function.\n",
    "\n",
    "We can plot this function to see what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAbQUlEQVR4nO3deZgU9bXG8e8Jxh0VBXFFIDFuUaN34jVGvRrU4BJJ4gKKigsSF1RURBBFXK5xV4zbJYjiFlBAxS2uGCIRdCCoKBpAUUEU3NBERZZz//jVJOM4Az0zXfXr7no/zzPP9HTXTL1UN2dqTlWfMndHRETy43uxA4iISLZU+EVEckaFX0QkZ1T4RURyRoVfRCRnVokdoBCtW7f29u3bx44hIlJWpkyZ8pG7t6l7f1kU/vbt21NdXR07hohIWTGzd+q7X60eEZGcUeEXEckZFX4RkZxR4RcRyRkVfhGRnEmt8JvZcDNbYGbTa923vpk9ZWYzk8+t0lq/iIjUL809/juAznXu6w884+5bAs8kX4uISIZSK/zuPgH4pM7dXYARye0RwK/TWr+ISDlzh0mT0vnZWff427r7/OT2B0DbhhY0s15mVm1m1QsXLswmnYhIibjxRjjlFFi6tPg/O9rBXQ9XgGnwKjDuPtTdq9y9qk2b77zjWESkYk2aBJdcAqNHwyopzFfIuvB/aGYbAySfF2S8fhGRkrZwIRx+OAwbBh07prOOrAv/OKBHcrsH8FDG6xcRKVnLlkH37nDkkXDwwemtJ83TOf8EvABsZWZzzewE4HJgXzObCeyTfC0iIsDFF8OSJXDppemuJ7XpnO5+RAMPdUprnSIi5erPf4bbboPq6nT6+rWVxVhmEZFK9s47cOyxcP/9sNFG6a9PIxtERCJavBgOOwz69oU99shmnSr8IiIRnXUWbLYZnH12dutUq0dEJJJ774Unnwx9fbPs1qvCLyISwWuvwRlnwNNPw7rrZrtutXpERDL2xRdwyCFw1VWw447Zr1+FX0QkQ+5wwgmw557hTJ4Y1OoREcnQDTfA7NkwcWK8DCr8IiIZ+dvf4LLLwhC21VePl0OtHhGRDCxYAF27hnfndugQN4sKv4hIypYtC4PXjj4aDjoodhoVfhGR1A0eDMuXhyFspUA9fhGRFD32GNx+O0yZkv7wtUKVSAwRkcozZw4cdxyMGQNtG7zQbPbU6hERSUHN8LVzz4Xdd4+d5ttU+EVEUtCnD2yxBZx5Zuwk36VWj4hIkd11FzzzTPbD1wqlwi8iUkSvvhpGLT/7LKyzTuw09VOrR0SkSD7/HA49FK69FrbfPnaahqnwi4gUgTscfzzsvXd4o1YpU6tHRKQIrr8+nL55992xk6ycCr+ISDNNnAiXXw6TJ8cdvlYotXpERJphwQLo1g2GD4f27WOnKYwKv4hIEy1bBkccAT16wIEHxk5TOBV+EZEmGjQonKd/0UWxkzSOevwiIk3wyCNw551h+FqLFrHTNI4Kv4hII731Vjh188EHYcMNY6dpPLV6REQa4euvw/C1gQNht91ip2kaFX4RkUY4/XT44Q/D53KlVo+ISIFGjIAJE+Cll0pz+FqhVPhFRArwyivQty889xy0bBk7TfNEafWY2Zlm9pqZTTezP5lZGbzXTUTyatEiOOSQMJZhu+1ip2m+zAu/mW0KnA5UufuPgRZAt6xziIgUomb42r77QvfusdMUR6xWzyrAGma2BFgTeD9SDhGRFbr2WnjvPbj33thJiifzwu/u88zsauBd4CvgSXd/su5yZtYL6AXQrl27bEOKiAB//StceSW8+CKstlrsNMUTo9XTCugCdAA2AdYys6PqLufuQ929yt2r2rRpk3VMEcm5Dz4Ic3juuCNcO7eSxDi4uw/wtrsvdPclwFigTN8GISKVaOnSMHHzhBNg//1jpym+GIX/XWBXM1vTzAzoBMyIkENEpF7nnx9aO4MGxU6Sjhg9/slmNhqYCiwF/g4MzTqHiEh9xo0LB3KnTi2/4WuFinJWj7tfCFwYY90iIg2ZPRt69gzFv3Xr2GnSo1k9IiLAV1/BoYfCBRfArrvGTpMuFX4REeC002DrraF379hJ0qdZPSKSe7ffHi6YXu7D1wqlwi8iuTZtGvTrB3/5C6y9duw02VCrR0Ry67PPQl//hhtg221jp8mOCr+I5JI7HHccdO4c3qGbJ2r1iEguXX01vP8+jBwZO0n2VPhFJHcmTIBrrqm84WuFUqtHRHJl/vzQ2hkxAvI6+FeFX0Ryo2b42oknwi9/GTtNPCr8IpIb550Ha6wR3p2bZ+rxi0guPPggjBoFU6ZU7vC1Qqnwi0jFmzULevWCRx6p7OFrhVKrR0QqWs3wtQsvhF12iZ2mNKjwi0hFO/XU8K7cU06JnaR0qNUjIhXrtttg0qRwvn4ehq8VSoVfRCrS3/8O/fuHN2vlZfhaodTqEZGKUzN87cYbYZttYqcpPSr8IlJRli+HHj3gwAOha9fYaUqTWj0iUlGuugoWLID774+dpHSp8ItIxRg/Hq6/PhzMXXXV2GlKl1o9IlIR3n8funeHO++EzTePnaa0qfCLSNlbsiT0808+GfbdN3aa0qfCLyJlb8AAaNkSBg6MnaQ8qMcvImVt7FgYPToMX/uedmULosIvImVr5kw46SR49FHYYIPYacqHfj+KSFn68ks45BC46CL46U9jpykvKvwiUnbcw9C1HXYIe/zSOGr1iEjZGTYMqqth8mQNX2sKFX4RKStTpoRLKD7/PKy1Vuw05SlKq8fM1jOz0Wb2hpnNMLOfxcghIuXlk0/gsMPglltgq61ipylfsfb4hwB/dvdDzWxVYM1IOUSkTCxfDsccA126hMmb0nSZF34zWxfYEzgWwN2/Ab7JOoeIlJfLL4dPP4Urr4ydpPzFaPV0ABYCt5vZ381smJl9p1NnZr3MrNrMqhcuXJh9ShEpGc8+C3/4A9x3H3z/+7HTlL8YhX8VYGfgFnffCfgX0L/uQu4+1N2r3L2qTZs2WWcUkRIxb14Yvnb33bDpprHTVIYYhX8uMNfdJydfjyb8IhAR+Zaa4Wu9e0OnTrHTVI7MC7+7fwC8Z2Y1x+Q7Aa9nnUNESt+558J664UhbFI8sc7qOQ24Jzmj5y3guEg5RKREjR4NDzyg4WtpKKjwJ+fZHwXsAWwMfAVMBx4F7nb3RY1ZqbtPA6oaF1VE8uLNN8Ns/ccfh/XXj52m8qz096iZPQ70BJ4AOhMK/7bA+cDqwENmdnCaIUUkP/71rzB87dJLoUq7h6koZI//aHf/qM59/wSmJh/XmFnroicTkdxxD0PXdt4ZevWKnaZyrXSPv6bom9kFZvatK1maWa/ay4iINMf//R+8/DLcequGr6WpMYdMTgP+bGZ717pPA1FFpCiqq2HQoHBQd00NcUlVYwr/PGB/4HIzOye5T7+TRaTZPv44zN+59Vb40Y9ip6l8jTpJyt3fBf4H2NbM7gfWSCWViOTG8uVw9NHhgO5vfxs7TT40pvBXA7j71+5+HPAcsGoaoUQkPy67DL74Igxhk2wUXPjd/cQ6X9/k7h2LH0lE8uLpp+Hmm2HUKA1fy1Ih5/E/bGa/MrPvPC1m1tHMLjaz49OJJyKVau7c0OK55x7YZJPYafKlkPP4TwTOAq43s08II5VXJ4xXngXc6O4PpRdRRCrNN9+EK2mdfjrsvffKl5fiWmnhT4aq9QP6mVl7/jOy4R/u/mWq6USkIp1zDrRuHYawSfYaNaTN3ecAc1JJIiK5MGoUPPywhq/FtNLCb2ZfAF7rLgc+AsYD57r7xyllE5EK88YbYbb+E09Aq1ax0+RXISMbWrr7OrU+1iVM1nwNuDX1hCJSEf75z3Cu/mWXhVk8Ek+T/tBy90/d/TrgB0XOIyIVyB1+9zv46U+hZ8/YaaTJF2JJTu+MdSEXESkjt9wC06fDCy9o+FopKKTHX9+bqFsBXQnXyxURadCLL8LgwTBxooavlYpC9th/VedrBz4Ghrj7o8WPJCKV4uOP4fDDw7jlLbeMnUZqFHIef4PXwzWzd929XXEjiUglWL4cjjoqvFHrN7+JnUZqa+5ZtOrWiUi9Lr00XEbx97+PnUTqau7BWV/5IiKSN08+Gdo71dWwik4BKTmFHNw9q6GHgLWLG0dEyt1778Exx8DIkbDxxrHTSH0K+V3ccgWPDSlWEBEpfzXD1848E/baK3YaaUghB3cvyiKIiJS/vn2hbVvo1y92ElmRgg/umtmPzOwZM5uefL2DmZ2fXjQRKScjR8Kjj8KIEXqTVqlrzFk9fwQGAEsA3P0VoFsaoUSkvMyYAaedBmPGwHrrxU4jK9OYwr+mu79Y576lxQwjIuWnZvjaFVfAT34SO40UojGF/yMz+wHJKZxmdigwP5VUIlIW3KFXL/jZz+B4XYC1bDTmDNtTgaHA1mY2D3gb6J5KKhEpCzffHNo8f/tb7CTSGAUXfnd/C9jHzNYi/KXwJaHH/05K2USkhE2eDBddFCZurrFG7DTSGCtt9ZjZOmY2wMxuNLN9CQW/B+FC64enHVBESs9HH4Xha3/8I/xAV+UoO4Xs8d8FfAq8AJwIDCS8a/c37j6tqSs2sxZANTDP3Q9q6s8RkWwtWwbdu0O3btClS+w00hSFFP6O7r49gJkNIxzQbefuXzdz3WcAM4B1mvlzRCRDl1wCixfD//5v7CTSVIWc1bOk5oa7LwPmNrfom9lmwIHAsOb8HBHJ1hNPhPbOyJEavlbOCnnqdjSzz5PbBqyRfG2Au3tT9tivB/qxgjlAZtYL6AXQrp1G/ovE9u670KMH3HcfbLRR7DTSHCvd43f3Fu6+TvLR0t1XqXW70UXfzA4CFrj7lJWsd6i7V7l7VZs2bRq7GhEposWLw/C1vn1hzz1jp5Hmau6FWJri58DBZjYHGAn8wszujpBDRAp09tmw6abhs5S/zLt07j6AMPMHM9sL6OvuR2WdQ0QKc++9obdfXa3ha5VCh2dEpEGvvw5nnAFPPw3rrhs7jRRL1MLv7s8Bz8XMICL1++KLMHztqqtgxx1jp5FiitHjF5ES5w4nngi77w7HHhs7jRSbWj0i8h033gj/+IeGr1UqFX4R+ZYXXoBLLw2fV189dhpJg1o9IvJvCxdC164wbBh07Bg7jaRFhV9EgDB87cgjwwC2X/0qdhpJkwq/iABhtv6yZWEIm1Q29fhFhMcfh+HDYcoUDV/LAz3FIjn3zjvhlM3Ro6Ft29hpJAtq9Yjk2OLFcOih0K8f7LFH7DSSFRV+kRw780xo1w7OOit2EsmSWj0iOXXPPWEGz0svafha3qjwi+TQ9OnQpw8884yGr+WRWj0iOfP552H42jXXwA47xE4jMajwi+SIO/TsCXvtBcccEzuNxKJWj0iODBkCs2fDxImxk0hMKvwiOTFxIvz+9zBpkoav5Z1aPSI5sGABdOsW3p3boUPsNBKbCr9IhasZvnbMMXDggbHTSClQ4RepcBdeGA7qXnxx7CRSKtTjF6lgjz4KI0aE4WstWsROI6VChV+kQs2ZA8cfD2PHwoYbxk4jpUStHpEK9PXXYfha//7w85/HTiOlRoVfpAL16RPO3unTJ3YSKUVq9YhUmLvugvHjNXxNGqbCL1JBXn01jFh+9llYZ53YaaRUqdUjUiEWLQrD1667DrbfPnYaKWUq/CIVwD2cwdOpExx1VOw0UurU6hGpANddB+++C/feGzuJlAMVfpEy9/zzcMUVMHkyrLZa7DRSDtTqESljH34Yhq/dfju0bx87jZSLzAu/mW1uZuPN7HUze83Mzsg6g0glWLoUjjgCjjsODjggdhopJzFaPUuBs919qpm1BKaY2VPu/nqELCJla9CgMH9n8ODYSaTcZF743X0+MD+5/YWZzQA2BVT4RQr08MNw990aviZNE7XHb2btgZ2AyfU81svMqs2seuHChVlHEylZb78drps7ciS0aRM7jZSjaIXfzNYGxgB93P3zuo+7+1B3r3L3qjZ6dYsA/xm+dt55sNtusdNIuYpS+M3s+4Sif4+7j42RQaQcnX46/PCH4bNIU2Xe4zczA24DZrj7tVmvX6Rc3XEHTJig4WvSfDH2+H8OHA38wsymJR86GU1kBV5+Gc45B8aMgZYtY6eRchfjrJ7nAe2viBRo0aLQ1x8yBLbbLnYaqQR6565ICXOHY4+F/faDI4+MnUYqhWb1iJSwa66B998Pp26KFIsKv0iJmjABrr5aw9ek+NTqESlBH3wQ5vDccQdssUXsNFJpVPhFSszSpWHiZs+e0Llz7DRSiVT4RUrM+eeH1s6gQbGTSKVSj1+khDz0ULiK1tSpGr4m6VHhFykRs2fDiSfCuHHQunXsNFLJ1OoRKQFffQWHHBLaO7vuGjuNVDoVfpES0Ls3bLMNnHpq7CSSB2r1iEQ2fDi88AK8+KKGr0k2VPhFIpo2Dc49N7xZa+21Y6eRvFCrRySSzz4Lw9f+8IfQ5hHJigq/SAQ1w9f23z+8WUskS2r1iERw1VVhLMN998VOInmkwi+Ssb/8Ba69NhzMXXXV2Gkkj9TqEcnQ/Plhrv6dd0K7drHTSF6p8ItkpGb4Wq9e4cIqIrGo8Itk5LzzYM014YILYieRvFOPXyQDDzwQDuROmQLf0+6WRKbCL5KymTPhd7+DRx6BDTaInUZErR6RVH35ZXiT1uDBsMsusdOIBCr8IilxD0PXfvxjOPnk2GlE/kOtHpGU3HZbOFdfw9ek1Kjwi6Rg6lQYMAD++ldYa63YaUS+Ta0ekSL79FM47DC46SbYeuvYaUS+S4VfpIiWL4cePeCgg+Dww2OnEamfWj0iRXTllfDRRzB6dOwkIg1T4RcpkvHjYcgQeOklDV+T0qZWj0gRvP8+dO8ehq9ttlnsNCIrpsIv0kxLlkDXruFc/X33jZ1GZOWiFH4z62xmb5rZLDPrHyODSLEMGAAtW8LAgbGTiBQm8x6/mbUAbgL2BeYCL5nZOHd/PessIs01diyMGQPV1Rq+JuUjxsHdXYBZ7v4WgJmNBLoARS/8Tz0FDz5Y7J8qEriHs3cee0zD16S8xCj8mwLv1fp6LvDfdRcys15AL4B2TbxUUevWsO22TfpWkYJ06wZVVbFTiDROyZ7O6e5DgaEAVVVV3pSfsdNO4UNERP4jRldyHrB5ra83S+4TEZEMxCj8LwFbmlkHM1sV6AaMi5BDRCSXMm/1uPtSM+sNPAG0AIa7+2tZ5xARyasoPX53fwx4LMa6RUTyTmcei4jkjAq/iEjOqPCLiOSMCr+ISM6Ye5PeG5UpM1sIvNPEb28NfFTEOMWiXI2jXI2jXI1Tqbm2cPc2de8si8LfHGZW7e4l96Z65Woc5Woc5WqcvOVSq0dEJGdU+EVEciYPhX9o7AANUK7GUa7GUa7GyVWuiu/xi4jIt+Vhj19ERGpR4RcRyZmKKPxmdpiZvWZmy82sqs5jA5KLur9pZr9s4Ps7mNnkZLlRybjoYmccZWbTko85ZjatgeXmmNmryXLVxc5Rz/oGm9m8WtkOaGC5zsk2nGVm/TPIdZWZvWFmr5jZA2a2XgPLZbK9VvbvN7PVkud4VvJaap9Wllrr3NzMxpvZ68nr/4x6ltnLzBbVen4HpZ0rWe8KnxcLbki21ytmtnMGmbaqtR2mmdnnZtanzjKZbC8zG25mC8xseq371jezp8xsZvK5VQPf2yNZZqaZ9WhSAHcv+w9gG2Ar4Dmgqtb92wIvA6sBHYDZQIt6vv8+oFty+1bg5JTzXgMMauCxOUDrDLfdYKDvSpZpkWy7jsCqyTbdNuVc+wGrJLevAK6Itb0K+fcDpwC3Jre7AaMyeO42BnZObrcE/lFPrr2AR7J6PRX6vAAHAI8DBuwKTM44XwvgA8IbnDLfXsCewM7A9Fr3XQn0T273r+81D6wPvJV8bpXcbtXY9VfEHr+7z3D3N+t5qAsw0t0Xu/vbwCzCxd7/zcwM+AUwOrlrBPDrtLIm6zsc+FNa60jBLsAsd3/L3b8BRhK2bWrc/Ul3X5p8OYlwpbZYCvn3dyG8diC8ljolz3Vq3H2+u09Nbn8BzCBc07ocdAHu9GASsJ6ZbZzh+jsBs929qRMBmsXdJwCf1Lm79muooTr0S+Apd//E3T8FngI6N3b9FVH4V6C+C7vX/Y+xAfBZrSJT3zLFtAfwobvPbOBxB540synJBeez0Dv5c3t4A39eFrId03Q8Ye+wPllsr0L+/f9eJnktLSK8tjKRtJZ2AibX8/DPzOxlM3vczLbLKNLKnpfYr6luNLzzFWN7AbR19/nJ7Q+AtvUsU5TtVrIXW6/LzJ4GNqrnoYHu/lDWeepTYMYjWPHe/u7uPs/MNgSeMrM3kr2DVHIBtwCXEP6jXkJoQx3fnPUVI1fN9jKzgcBS4J4GfkzRt1e5MbO1gTFAH3f/vM7DUwntjH8mx28eBLbMIFbJPi/JMbyDgQH1PBxre32Lu7uZpXaufdkUfnffpwnfVsiF3T8m/Jm5SrKn1uSLv68so5mtAvwW+K8V/Ix5yecFZvYAoc3QrP8whW47M/sj8Eg9DxWyHYuey8yOBQ4COnnS4KznZxR9e9WjkH9/zTJzk+d5XcJrK1Vm9n1C0b/H3cfWfbz2LwJ3f8zMbjaz1u6e6kCyAp6XVF5TBdofmOruH9Z9INb2SnxoZhu7+/yk7bWgnmXmEY5D1NiMcGyzUSq91TMO6JaccdGB8Jv7xdoLJAVlPHBoclcPIK2/IPYB3nD3ufU9aGZrmVnLmtuEA5zT61u2WOr0VX/TwPpeAra0cPbTqoQ/k8elnKsz0A842N2/bGCZrLZXIf/+cYTXDoTX0rMN/bIqluQYwm3ADHe/toFlNqo51mBmuxD+z6f6C6nA52UccExyds+uwKJabY60NfhXd4ztVUvt11BDdegJYD8za5W0ZfdL7muctI9eZ/FBKFhzgcXAh8ATtR4bSDgj401g/1r3PwZsktzuSPiFMAu4H1gtpZx3ACfVuW8T4LFaOV5OPl4jtDzS3nZ3Aa8CryQvvI3r5kq+PoBw1sjsjHLNIvQypyUft9bNleX2qu/fD1xM+MUEsHry2pmVvJY6ZrCNdie06F6ptZ0OAE6qeZ0BvZNt8zLhIPluGeSq93mpk8uAm5Lt+Sq1zsZLOdtahEK+bq37Mt9ehF8884ElSe06gXBM6BlgJvA0sH6ybBUwrNb3Hp+8zmYBxzVl/RrZICKSM5Xe6hERkTpU+EVEckaFX0QkZ1T4RURyRoVfRCRnVPhFRHJGhV9EJGdU+EWawMxOqjWz/W0zGx87k0ih9AYukWZIZuU8C1zp7g/HziNSCO3xizTPEMJcHhV9KRtlM51TpNQk00O3IMx3ESkbavWINIGZ/RfhKkl7eLgSkkjZUKtHpGl6E657Oj45wDssdiCRQmmPX0QkZ7THLyKSMyr8IiI5o8IvIpIzKvwiIjmjwi8ikjMq/CIiOaPCLyKSM/8PfOVuO6qo0pMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = range(-10, 10+1, 1)\n",
    "y = [activation_function(z) for z in x]\n",
    "plt.plot(x, y, c='blue', lw='1')\n",
    "plt.ylabel('ReLU(z)')\n",
    "plt.xlabel('z')\n",
    "plt.show()\n",
    "plt.clf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note it's a combination of linear ($f(x) = x$) and constant ($f(x) = 0$) functions. A negative signal is clamped to zero but otherwise unperturbed.\n",
    "\n",
    "Now we've determined how a neuron will accept an input and produce an output. These two functions are all we need to implement our feedforward function.\n",
    "\n",
    "### Feedforward Function\n",
    "\n",
    "For our feedforward function, we only needed the summation and activation functions as previously defined. Data is fed into the network starting at the first layer, it is processed by the neurons in the input layer, neural outputs are calculated as a function of the data, the neuron's weighted sum of the data, and this value transferred to the activation function. Then, the next layer takes those outputs as inputs and the process cycles until you reach the output layer.\n",
    "\n",
    "Since we don't need to store outputs within the neurons (as would be the case for backpropagation), the feedforward function need only return a value and not mutate the network in any way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feed_forward(network, example):\n",
    "    \"\"\"Feedforward method. Feeds data forward through network.\n",
    "    \n",
    "    Parameters:\n",
    "        network : the neural network.\n",
    "        example : an example of data to feed forward.\n",
    "    \n",
    "    Returns:\n",
    "        The output of the forward pass.\n",
    "    \"\"\"\n",
    "    layer_input, layer_output = example, []\n",
    "    for layer in network:\n",
    "        for neuron in layer:\n",
    "            # sum the weight with inputs\n",
    "            summ = summing_function(neuron, layer_input)\n",
    "            # activate the sum, append output to outputs\n",
    "            layer_output.append(activation_function(summ))\n",
    "        # inputs become outputs of previous layer\n",
    "        layer_input, layer_output = layer_output, []\n",
    "    return layer_input # return the final output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've done all of this because our error function demands the use of data being fed into the network to calculate the network error.\n",
    "\n",
    "### Mean Squared Error\n",
    "\n",
    "Now that we can feedforward data through the network, we can implement our error function, the mean squared error. The mean squared error function will feed all training data into the network and find the average error between them. The error is calculated as a function of what is known about the data (the classification is already known) versus what the network outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sse(actual, target):\n",
    "    \"\"\"Sum of Squared Error.\n",
    "    \n",
    "    Parameters:\n",
    "        actual : network output.\n",
    "        target : example target output.\n",
    "    \n",
    "    Returns:\n",
    "        The sum of squared error of the network for an example.\n",
    "    \"\"\"\n",
    "    summ = 0.00\n",
    "    for i in range(len(actual)):\n",
    "        summ += (actual[i] - target[i])**2\n",
    "    return summ\n",
    "\n",
    "def mse(network):\n",
    "    \"\"\"Mean Squared Error.\n",
    "    \n",
    "    Parameters:\n",
    "        network : the neural network to test.\n",
    "    \"\"\"\n",
    "    training = TRAIN\n",
    "    summ = 0.00\n",
    "    # for each training example\n",
    "    for example in training:\n",
    "        # populate a target vector\n",
    "        target = [0 for _ in range(CLASSES)]\n",
    "        # denote correct classification\n",
    "        target[int(example[-1])] = 1\n",
    "        # get actual output by feeding example through network\n",
    "        actual = feed_forward(network, example)\n",
    "        # sum up the sum of squared error\n",
    "        summ += sse(actual, target)\n",
    "    # MSE is just sum(sse)/number of examples\n",
    "    return summ / len(training)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the mean squared error is simply the mean of the sum of squared error per training example. This gives us the network error considering every training example.\n",
    "\n",
    "## The Particle Swarm Part\n",
    "\n",
    "So far, we've only been fooling around with the neural network side of things when the brunt of this tutorial will instead focus upon the PSO side of things. At this point we have our fitness function and can begin to populate a swarm with particles.\n",
    "\n",
    "As previously mentioned, we're using a PSO algorithm to train the network. This is a misnomer, as the use of a metaheuristic such as PSO instead tries to find weights and doesn't update weights as an informed decision. It is, however, an informed search. We look at a swarm of possible weights, see which works the best, and then update the swarm based on the swarm best. To compare, a random search wouldn't be informed whereas a PSO algorithm makes decisions which ensure we get a better fitness over time.\n",
    "\n",
    "### Representing Particles\n",
    "\n",
    "A particle for our algorithm is some encoding of network weights. Given a network, we need a way to represent it as a single array or list. Thankfully it is relatively straightforward: we take our network, take each weight in order, and represent it as a point in `n`-space, where `n` is the number of weights and biases. These points are our particles. For the Wine data set, our particles are in `105`-space.\n",
    "\n",
    "First we'll make a class to containerize a particle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Particle:\n",
    "    \"\"\"Particle class.\n",
    "    Containzerizes a position, velocity.\n",
    "    \n",
    "    Attributes:\n",
    "        pos : the position in n-space.\n",
    "        best_pos : the best position this particle has had.\n",
    "        fit : the fitness of the particle.\n",
    "        best_fit : the best fitness this particle has had.\n",
    "        vel : the velocity in n-space.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, pos, vel):\n",
    "        \"\"\"Particle constructor.\"\"\"\n",
    "        # initialize position and velocity as params\n",
    "        self.pos, self.vel = pos, vel\n",
    "        # find fitness at instantiation\n",
    "        network = initialize_network(self.pos)\n",
    "        self.fit = mse(network)\n",
    "        # best so far is just initial\n",
    "        self.best_pos, self.best_fit = self.pos, self.fit\n",
    "    \n",
    "    def set_pos(self, pos):\n",
    "        \"\"\"Position mutator method.\"\"\"\n",
    "        self.pos = pos\n",
    "        if not any(p < -7.00 for p in pos)\\\n",
    "        and not any(p > 7.00 for p in pos):\n",
    "            # get fitness of new position\n",
    "            network = initialize_network(self.pos)\n",
    "            fitness = mse(network)\n",
    "            # if better\n",
    "            if fitness < self.best_fit:\n",
    "                self.fit = fitness\n",
    "                # update best fitness\n",
    "                self.best_fit = self.fit\n",
    "                # update best position\n",
    "                self.best_pos = self.pos\n",
    "    \n",
    "    def set_vel(self, vel):\n",
    "        \"\"\"Velocity mutator method.\"\"\"\n",
    "        self.vel = vel\n",
    "    \n",
    "    def get_pos(self):\n",
    "        \"\"\"Position accessor method.\"\"\"\n",
    "        return self.pos\n",
    "    \n",
    "    def get_vel(self):\n",
    "        \"\"\"Velocity accessor method.\"\"\"\n",
    "        return self.vel\n",
    "    \n",
    "    def get_best_pos(self):\n",
    "        \"\"\"Best position accessor method.\"\"\"\n",
    "        return self.best_pos\n",
    "\n",
    "    def get_fit(self):\n",
    "        \"\"\"Fitness accessor method.\"\"\"\n",
    "        return self.fit\n",
    "    \n",
    "    def get_best_fit(self):\n",
    "        \"\"\"Best fitness accessor method.\"\"\"\n",
    "        return self.best_fit\n",
    "    \n",
    "    def __lt__(self, other):\n",
    "        return self.fit < other.fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we now know how particles are represented, we can make a swarm of them.\n",
    "\n",
    "### Swarm Initialization\n",
    "\n",
    "First, we need to initialize an initial swarm. Since we don't know our search space bounds, we'll use an initialization strategy placing each particle within `[-1.00, 1.00]`, our swarm consists of particles in `105`-space where each axis is randomized within this interval. Recall a particle is just an encoding of some solution to the problem: our solutions are in `105` dimensions so particles are also in `105` dimensions. Additionally, we'll start with a swarm size of `100`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "SWARM_SIZE = 100 # how large our swarm is\n",
    "\n",
    "def initialize_swarm(size, dim):\n",
    "    \"\"\"Swarm initialization function.\n",
    "    \n",
    "    Parameters:\n",
    "        size : the size of our swarm.\n",
    "        dim : the dimensionality of the problem.\n",
    "    \n",
    "    Returns:\n",
    "        A random swarm of that many Particles.\n",
    "    \"\"\"\n",
    "    swarm = [] # swarm stored as list\n",
    "    for _ in range(size): # for the size of the swarm\n",
    "        # position is random in every dimension\n",
    "        position = [random.uniform(-1.00, 1.00) for _ in range(dim)]\n",
    "        # velocity is initially zero in every dimension\n",
    "        velocity = [0 for _ in range(dim)]\n",
    "        # init a particle\n",
    "        particle = Particle(position, velocity)\n",
    "        swarm.append(particle) # add to swarm\n",
    "    return swarm\n",
    "\n",
    "SWARM = initialize_swarm(SWARM_SIZE, DIMENSIONS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't need to sort our population by fitness but we'll still make a function to do this and print the five best particle fitnesses to verify our algorithm works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Particle 1 = 0.987354856821158\n",
      "Particle 2 = 1.07611596069502\n",
      "Particle 3 = 1.1569260312369962\n",
      "Particle 4 = 1.2170419249411921\n",
      "Particle 5 = 1.241404249558855\n"
     ]
    }
   ],
   "source": [
    "def print_five_best():\n",
    "    SWARM.sort()\n",
    "    for i in range(5):\n",
    "        print(f'Particle {i+1} = {SWARM[i].get_fit()}')\n",
    "\n",
    "print_five_best()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use this function a couple more times to verify that solutions improve over time.\n",
    "\n",
    "### Moving Particles\n",
    "\n",
    "Particles move within the search space based on the PSO parameters: the inertial weight and social and cognitive coefficients. These should be empirically found for each problem and shouldn't be known a priori but there are a few shortcuts that can be used: some parameter choices have been empirically found to work well if not great for most problems.\n",
    "\n",
    "Here we'll use some parameters I've found to work well for this specific network and data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = 0.3 # intertial weight\n",
    "C_1 = 1.5 # cognitive coefficient\n",
    "C_2 = 1.2 # social coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To actually move the particles, we need a way to find the swarm best position:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_swarm_best(swarm):\n",
    "    \"\"\"Finds the swarm best fitness and position.\n",
    "    \n",
    "    Parameters:\n",
    "        swarm : the swarm to search.\n",
    "        \n",
    "    Returns:\n",
    "        The swarm best fitness and swarm best position.\n",
    "    \"\"\"\n",
    "    # initially assume the first is the best\n",
    "    best_fit = swarm[0].get_fit()\n",
    "    best_pos = swarm[0].get_pos()\n",
    "    for particle in swarm: # for each particle\n",
    "        # if better fitness found\n",
    "        if particle.get_fit() < best_fit:\n",
    "            # update best fitness and position\n",
    "            best_fit = particle.get_fit()\n",
    "            best_pos = particle.get_pos()\n",
    "    return best_fit, best_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use this function several times when we update our positions. And about updating the positions, it is simply a function of the swarm best position, each particle's best position, their current positions and velocities, and the PSO parameters as defined previously. The math is below:\n",
    "\n",
    "$$p_i(t+1) = p_i(t) + v_i(t+1)$$\n",
    "\n",
    "Or, for a dimension $i$, the new position in that dimension is a function of its prior position and new velocity. Here $(t+1)$ means the current iterative step and $t$ is the prior. The new velocity is calculated as:\n",
    "\n",
    "$$v_i(t+1) = \\omega v_i(t) + c_1r_1(y(t) - x_i(t)) + c_2r_2(\\hat{y}(t) - x_i(t))$$\n",
    "\n",
    "Where $\\omega$ is the inertial weight, $c_1$ and $c_2$ are the cognitive and social coefficients respectively, $r_1$ and $r_2$ are stochastic multiplicands, and $y$ and $\\hat{y}$ is the particle best position and swarm best position respectively. It's useful to think of it split into three separate terms: a weight term ($w$), a cognitive term ($c$) and a social term ($s$):\n",
    "\n",
    "$$\\begin{align}\n",
    "w &= \\omega v_i(t) \\nonumber \\\\\n",
    "c &= c_1r_1 (y(t) - x_i(t)) \\nonumber \\\\\n",
    "s &= c_2r_2 (\\hat{y}(t) - x_i(t)) \\nonumber \n",
    "\\end{align}$$\n",
    "\n",
    "Since each term in order deals with the inertial weight, cognitive aspect, and social aspect of the velocity update. The new velocity is simply:\n",
    "\n",
    "$$v_i(t+1) = w + c + s$$\n",
    "\n",
    "We'll implement particle movement below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def move_particles(swarm, dim, ic, cc, sc):\n",
    "    \"\"\"Particle movement function.\n",
    "    \n",
    "    Parameters:\n",
    "        swarm : the swarm to move.\n",
    "        dim : dimensionality of the network.\n",
    "        ic : inertial coefficient.\n",
    "        cc : cognitive coefficient.\n",
    "        sc : social coefficient.\n",
    "    \"\"\"\n",
    "    # get swarm bests\n",
    "    swarm_best = get_swarm_best(swarm)\n",
    "    for particle in swarm: # for each particle\n",
    "        # new position and velocity is initially zero\n",
    "        new_pos = [0 for _ in range(dim)]\n",
    "        new_vel = [0 for _ in range(dim)]\n",
    "        for d in range(dim): # for each axis\n",
    "            # the social and cognitive coefficients \n",
    "            # take a stochastic multiplicand\n",
    "            r_1 = random.uniform(0.00, 1.00)\n",
    "            r_2 = random.uniform(0.00, 1.00)\n",
    "            # this is split for readability but the update is based\n",
    "            # on an addition of a weight, cognitive, and social term\n",
    "            weight = ic * particle.get_vel()[d]\n",
    "            cognitive = cc * r_1\n",
    "            cognitive *= (particle.get_best_pos()[d] - particle.get_pos()[d])\n",
    "            social = sc * r_2\n",
    "            social *= (swarm_best[1][d] - particle.get_pos()[d])\n",
    "            # new velocity is simply weight + cognitive + social\n",
    "            new_vel[d] = weight + cognitive + social\n",
    "            # new position is just old position + velocity\n",
    "            new_pos[d] = particle.get_pos()[d] + new_vel[d]\n",
    "        # update particle with new position and velocity\n",
    "        particle.set_pos(new_pos)\n",
    "        particle.set_vel(new_vel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then move the particles five times to see how the best fitness improves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial fitness: 0.9874\n",
      "Iteration 1 fitness: 0.6415\n",
      "Iteration 2 fitness: 0.6178\n",
      "Iteration 3 fitness: 0.5488\n",
      "Iteration 4 fitness: 0.5227\n",
      "Iteration 5 fitness: 0.4667\n"
     ]
    }
   ],
   "source": [
    "SWARM.sort()\n",
    "print(f'Initial fitness: {SWARM[0].get_fit():.4f}')\n",
    "\n",
    "for i in range(5):\n",
    "    move_particles(SWARM, DIMENSIONS, W, C_1, C_2)\n",
    "    best = get_swarm_best(SWARM)\n",
    "    print(f'Iteration {i+1} fitness: {best[0]:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the fitness was initially around `0.99` and moving the particles five times made the most fit particle move around such that the best fitness is now roughly `0.47`, our anticipated behavior.\n",
    "\n",
    "## Measuring Performance\n",
    "\n",
    "The goal to using PSO to train our network is to optimize the weights such that the mean squared error function is minimized. We know classification accuracy is (roughly) inversely proportional to the mean squared error but we should examine this as well. We'll use the below functions to make a nice graph later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def performance_measure(particle, data):\n",
    "    \"\"\"Measures accuracy of the network using classification error.\n",
    "\n",
    "    Parameters:\n",
    "        particle : the particle to test.\n",
    "        data : a set of data examples.\n",
    "    Returns:\n",
    "        A percentage of correct classifications.\n",
    "    \"\"\"\n",
    "    network = initialize_network(particle)\n",
    "    correct, total = 0, 0\n",
    "    for example in data:\n",
    "        # check to see if the network output matches target output\n",
    "        if check_output(network, example) == float(example[-1]):\n",
    "            correct += 1\n",
    "        total += 1\n",
    "    return 100*(correct / total)\n",
    "\n",
    "def check_output(network, example):\n",
    "    \"\"\"Compares network output to actual output.\n",
    "\n",
    "    Parameters:\n",
    "        network : the neural network.\n",
    "        example : an example of data.\n",
    "    Returns:\n",
    "        The class the example belongs to (based on network guess).\n",
    "    \"\"\"\n",
    "    output = feed_forward(network, example)\n",
    "    return output.index(max(output))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is just another way for us to look at the network performance. We could even use the classification accuracy as our PSO algorithm fitness function but it is safer to use the mean squared error for this reason instead.\n",
    "\n",
    "## Main Driver\n",
    "\n",
    "In essence, we've finished our PSO algorithm for training a network but we should make a main driver to move particles over several iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPOCHS = 100\n",
    "\n",
    "MSE = [] # to store MSE over epochs\n",
    "TRP = [] # training accuracy over epochs\n",
    "TEP = [] # testing accuracy over epochs\n",
    "\n",
    "def pso(dim, epochs, swarm_size, ic, cc, sc):\n",
    "    \"\"\"Particle Network training function\n",
    "    Main driver for PSO algorithm\n",
    "\n",
    "    Parameters:\n",
    "        dim : dimensionality of the problem.\n",
    "        epochs : how many iterations.\n",
    "        swarm_size : how big a swarm is.\n",
    "        ic : inertial coefficient (omega).\n",
    "        cc : cognitive coefficient (c_1).\n",
    "        sc : social coefficient (c_2).\n",
    "    \n",
    "    Returns:\n",
    "        Idealized network weights.\n",
    "    \"\"\"\n",
    "    swarm = initialize_swarm(swarm_size, dim) # init swarm\n",
    "    for e in range(1, epochs+1):\n",
    "        # get swarm best fitness and position\n",
    "        swarm_best = get_swarm_best(swarm)\n",
    "        MSE.append(swarm_best[0]) # get error of network using swarm best\n",
    "        # get classification error of network for training and test\n",
    "        TRP.append(performance_measure(swarm_best[1], TRAIN))\n",
    "        TEP.append(performance_measure(swarm_best[1], TEST))\n",
    "        # reposition particles based on PSO params\n",
    "        move_particles(swarm, dim, ic, cc, sc)\n",
    "    return get_swarm_best(swarm)[0]\n",
    "\n",
    "NETWORK = pso(DIMENSIONS, EPOCHS, SWARM_SIZE, W, C_1, C_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have three global lists, MSE, TRP, and TEP to store the mean squared error, training performance, and testing performance per epoch of the algorithm. We can plot these:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEaCAYAAACmbNjHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3yV9fn4/9eVQRKSQNgrIGEJQSNIZIhaHBQE1FoXqBVXqRZFqX5c/dUQW7/V2rq1iOOj1IFU6wesOChCwcUSZQSREILsESAEwsi4fn/cJ+Fkn4z7nJyc6/l4nEfu8b7vcx2O5sp73O+3qCrGGGNMMAkLdADGGGNMbVnyMsYYE3QseRljjAk6lryMMcYEHUtexhhjgo4lL2OMMUHHkpcxxphKichrIrJHRNZ6HWstIvNFZKPnZyvPcRGRZ0UkU0RWi8iZrsYWbM95hYWFaUxMTKDDMMaYoJefn6+qWmUlRkTOAw4DM1X1NM+xvwD7VfUxEXkAaKWq94vIGOBOYAwwBHhGVYe4FXuEWzd2S0xMDEeOHAl0GMYYE/RE5Gh151V1sYh0L3f4MmCEZ/sNYBFwv+f4THVqRN+ISIKIdFLVnQ0ZcwlXmw1FZLSIbPBUIx+o5PwpIrLAU8VcJCKJbsZjjDGm3jp4JaRdQAfPdhdgq1e5bZ5jrnAteYlIOPACcDGQDEwQkeRyxf6Kk6lTgEeAP7sVjzHGmAoiRGSF12tSbS721LIC0vfkZrPhYCBTVbMARGQWTrUyw6tMMvA7z/ZC4P9cjMcYY0xZhaqaWstrdpc0B4pIJ2CP5/h2oKtXuUTPMVe42WzoSxXye+CXnu3LgXgRaVP+RiIyqeQvg8LCQleCNcYY45O5wETP9kRgjtfxGzyjDocCuW71d0Hgh8rfC/xMRFYBP8PJ0kXlC6nqDFVNVdXUiIigG2NijDFBSUTeAb4GThWRbSJyC/AYMFJENgIXefYB5gFZQCbwMvBbN2NzMxPUWIVU1R14al4iEgdcoaoHXYzJGGOMj1R1QhWnLqykrAKT3Y3oJDdrXsuB3iKSJCLNgPE41cpSItJWREpieBB4zcV4jDHGNBGuJS9VLQTuAD4F1gOzVXWdiDwiIpd6io0ANojIjzjDLR91Kx5jjDFNR9DNsBEbG6t1fUh5xgy49VYIC3RPnzHGNAIikq+qsYGOoy5C6tf41KmQnx/oKGpPVXnmm2c4eMy6A40xBkIsecXHQ15eoKOonZz8HL7b9R0RYREcLah2JhdjjAkZIZe8Dh0KdBSVU1U+3PAhr3z7CoXFJ59lW7d3HXM2zGHy4Mk0j2zOki1LAhilMcY0DiH10FSLFo2z5rUxZyNTPplC9sFsOsZ15Lllz/H8xc/TrWU3zjvlPM475TwAsg9m83Hmx5x7yrkBjtgYYwIrpJJXY2w2/Mf3/2Dqp1O5f/j93DX0LrYc3MKfFv+JCe9PIO9EHgM7DiQ6IrrMNQNfGki4hNO2eVs27t/IKS1PobC4kD1H9tA9oTs7D+8kKjyK1jGt2bBvA71a9+Jo4VH2H91Pt5bd2J63ndjIWBKiE1i/bz192/Ql70Qeh44fIrFFIlsPbaVlVEvimsXxY86P9G3bl4PHDpJfkE/n+M5syd1Cm5g2REdEk3Ugiz5t+pBzNIeCogI6xnVk84HNdIzrSHhYOD/l/kSv1r3Ym78XVaV9bHs2HdhEYnwiirIjbwc9WvVg95Hd9pnsM9lnqsNn+nDCh0SGRwboN1jghFzyakzNhvuP7ueez+7h84mfk9IhhR15O4iKiOKyvpfxwtgX+GLLFyAVr8vYm0FkWCS92/Rm/qb5DO82nKMFR1m/bz3ndDuH1btXE98snqRWSXyc+TEXJl3IgaMH2HxwM0MTh/Ltzm9p17wdXVt25cMNHzK2z1h2Hd7F9kPbOavLWSzbvozEFol0iO3Ax5kfM67POLbkbmH/0f0M7DiQr7Z+Re/WvWkR1YKF2QsZ3Ws0WQeyOHLiCKd3OJ0lW5bQv31/osKj+Hrb11zU4yI25GygWIvp17YfCzcvZFDnQagqq3atYkT3EfaZ7DPZZ6rjZwoPC/f/L69GIKSGyl93HVx8MVx/fQMHVUe/+/R35BfkM33cdACuff9abk+93ZoFjTF+EcxD5UOu5tVYmg0z92fyxvdvkPHbk5Psv33F2wGMyBhjgoeNNgyQB/7zAL8b+js6xDnruK3auYp/fP+PAEdljDHBIaRqXoEcbXj4xGH25e8DYO2etSzbvoyZl88sPR8dEU3rmNaBCc4YY4JMSCWv+HjYt8//73us8Bgpf0+hSIsQBBHh6dFP0zyyeWmZPm360K9dP/8HZ4wxQSjkklcgal7PfPMMZ3Q8gw+u+aDS80XFRSQ+lcjGOzcS1yzOz9EZY0zwCbnk5e8+r71H9vLEV0/w1S1fVVkmPCyczDsziW0WlIN+jDHG70JqwEYg+rymLZrGdadfR582faos89XWr9iSu8WPURljTHALuZqXP5PX+r3rmZ0xmx8m/1BtuZIn6JPbJfspMmOMCW4hl7z81Wyoqtzz2T08MPwB2jRvU23Z61Ku809QxhjTRFizoUv+35L/x778fdwx+I5qy23Yt4ErZ1/pn6CMMaaJCKnk1VDNhpv2b+L++feXPrdV3nsZ7/HSypeYM34OURFR1d6re0J3/nTBn+oflDHGhJCQS16HDkF9p3Oc+ulUvtj6BckvJDN9xXSKiotKz63csZLbP7qdOePn0Cm+U7X3yTuex/vr36dv2771C8gYY0KMq31eIjIaeAYIB15R1cfKne8GvAEkeMo8oKrz3IqnWTOIiIBjxyAmpm73WLh5IWv3rGX95PVsyNnA5HmTefzLx0tnx8g+mM3Ll7zMwE4Da7xXztEcsg9m1y0QY4wJYa7NKi8i4cCPwEhgG7AcmKCqGV5lZgCrVPXvIpIMzFPV7tXdtz6zygO0awfr1kH79rW/tliLSZ2RygPnPMDV/a8GnIEZa/asoaCoAIAWUS3o3aZ3jfdatXMV/dr1q7BWlzHG+EswzyrvZrPhYCBTVbNU9QQwC7isXBkFWni2WwI7XIwHqF+/15ur3yQ6Ipqrkq8qPSYipHRIYVDnQQzqPMinxAXw0sqXyNibUXNBY4wxFbjZbNgF2Oq1vw0YUq7MNOAzEbkTiAUucjEeoO7D5fML8vn9579n9pWz6/X+xVrMvvx9pWt4GWOMqb1AD9iYALyuqonAGOAfIlIhJhGZJCIrRGRFYWFhvd6wrsPl3894nzM6nMGwrsN48usnaf9EHdodga+3fs1vP/ptna41xhjjcLPmtR3o6rWf6Dnm7RZgNICqfi0i0UBbYI93IVWdAcwAp8+rPkHVtdkw+2A2Z3Q4A4ApQ6awbu86irWYsIq5tkoFRQUM7zacoYlDax+AMcaYUm7WvJYDvUUkSUSaAeOBueXK/ARcCCAi/YBoYK+LMdW52XB73nY6x3cGYMHmBTw/5vlaJS5VZeQ/RpKxN4PwsPDaB2CMMaaUa8lLVQuBO4BPgfXAbFVdJyKPiMilnmL3AL8Wke+Bd4Ab1a3hjx51bTbcnredLi26ADBr7Szum38fH2/82OfrRYT3r37f5i80xpgG4OpzXp5ntuaVO/aw13YGMNzNGMqra7PhjrwddIl3ktfrv3idrANZPq98/Nmmz/hP1n/4y8i/1P6NjTHGVBBSE/NCPZoNDzk1rwVZC/h257fcftbt7Dmyh4TohBqvPafbOZzS8pQ6RGuMMaYygR5t6Hd1aTYsKCpg/9H9dIjtQL92/RjZcySrdq7iiS+fqPa6Yi3m3s/u5dDxQ5za9tR6RG2MMcZbSNa8apu8dh7eSbvYdoSHhRMu4Zze/nTCw8I595Rzq71OEFI7p9ImpvolUYwxxtROyNW86pK8vPu7Jrw/gQ05GwB4aMFD5BfkV3pNfkE+8zbOY/xp44kMj6xXzMYYY8oKueTVokXt+7xK+rsAPp/4eemIwaSEpDIzynvbdXgXi7IX1SdUY4wxVQi55FWXmtf2vO10ie/Cpv2beHP1m6XHfz3o18RHxVd6TY9WPXji59X3iRljjKkbS14+2H7IeUC5sLjs1FSP/PcRnvz6yQrl9x/dz5BXhlCsxfUJ1RhjTBVCcsBGrZsN87aT3C6ZU9ueWmbU4NShU4mJrLgwWKvoVrz1y7dqNQOHMcY0RiIyFbgVZxWQNcBNQCeclULaACuBX3lWD/GbkPvtWpeh8jvydtClRReumH0FS7ctLT2uKEu2LKlQfu6GuSS2SKxvqMYYE1Ai0gWYAqSq6mk4iwaPBx4HnlLVXsABnHlq/Srkkld9+rymj51OSoeU0uOHjh/inxn/LFO2sLiQf/3wL6t1GWOaigggRkQigObATuAC4D3P+TeAX/g7qJD7DRsdDYWFcMLHCq6qsv3QduKj4sncn1mmmTCxRSIvjn2xTPmIsAje+MUbNAtv1pBhG2OMGyJKlpvyvCZ5n1TV7cBfcSZR3wnk4jQTHvTMXwvOWo1d/Bk0hGDyEqld0+Gh404H2aFjh3h7zdsVzt/z6T1k7s8s3R/z1hiyDmQ1SKzGGOOyQlVN9XrN8D4pIq2Ay4AkoDPOosGjAxBnBSGXvKB2TYcls8knt0/muTHPVTh/yamXlJmg98lRT9I9oXsDRWqMMQF1EbBZVfeqagHwL5zJ1BM8zYhQ+VqNrrPkVYOS2TXeXP0m//7x3xXOn9vtXJqFN2PDvg1MWzSNvm37Wn+XMaap+AkYKiLNRURw1l/MABYCV3rKTATm+DuwkPwtW5tZNkpm10hul0yPVj0qnJ+xcgZ//O8fadu8LcMShzVwpMYYEziquhRnYMa3OMPkw3BWtb8f+J2IZOIMl3/V37GF3HNeUPtmw85xnTm9/emVzlF4W+ptzPx+JoXFhYzqNaqBIzXGmMBS1TQgrdzhLGBwAMIpFZI1r1olL0/Na+irQ/l+1/cVzosIuw7vsqZCY4zxo5CtefnabLjj8A4uSLqAFb9egaKVlrn/nPsbMDpjjDE1CcnqQm2Gym8/tJ32se15e83bVrsyxphGIiR/G9e2zyshOoGl25fWXNgYY4xfuJq8RGS0iGwQkUwReaCS80+JyHee148ictDNeEr42mxYWFzIniN76Nu2L89e/Kz7gRljjPGJa8lLRMKBF4CLgWRggogke5dR1amqOkBVBwDP4TwA5zpfmw13H95Nm5g2LN6ymEcXP+p+YMYYY3zi5oCNwUCmqmYBiMgsnGlGMqooP4GKwzEbjKqy8/BOOsd39rnZsGQ2+dPan0b72PZuhWaMMaaW3Gw27AJs9dqvcvJGETkFZ+6sz6s4P6lk4sjCwsLKitRIUXo924vDJw77nLxKZpOPDI+kd5vedXpfY4wxDa+xDNgYD7ynqkWVnVTVGSUTR0ZE1K2yGCZh9GjVg6wDWT7PsFGygvKjix/l/Yz36/S+xhhjGp6bzYbbga5e+9VN3jgemOxiLAD0at2LzP2ZdItPqbHmdeTEEV5d9Sp3Dr6Tmwbe5HZoxhhjasHNmtdyoLeIJIlIM5wENbd8IRHpC7QCvnYxFuBk8qqp2bBYi7nh/24gpUMKNw64kekrpnPwmF8GQhpjjPGBa8nLs1DZHcCnwHpgtqquE5FHRORSr6LjgVmqWvn0FQ2oZ6ueZO7PrLHZ8A+f/4E9R/bw0riXEBH25e+zB5SNMaYRET/kjAYVGxurR44cqdO18zfN57EvH2POLxfQoQNUdpt3177LgwseZOmtS2kX266e0RpjTOMlIvmqGhvoOOoipKoTJc2GsbFw7BgUVTI85OVvX+apUU+VJq7sg9mMfXusnyM1xhhTnZBKXl1bdmX34d0cLzpGbGzFfq9iLWbFjhWc3fXs0mOd4jrx1Kin/BypMcaY6oRU8ooIi6Bby25sPrCZFi3gwIGy5zfmbKRVTKsyzYW5x3NpHtncz5EaY4ypTkglL3CaDjcd2MTll8PPfgZvvgnFxc655TuWc1bns8qUX5S9iJnfzwxApMYYY6oScut5lfR7PfccXHMN3HMPPP00zJsHy7YvY3CXsouDXt3/6gBFaowxpiohV/MqGS4PcM458PXXzkS9S5dWXvN6bdVrrNm9JhChGmOMqULIJa+SmleJsDDo3Bn2HTjB6t2rGdR5UJnybWLaEBMZ4+8wjTHGVCMkmw03HdhU5ljLlvDD/rUkJSQR1yyuzLnL+l7mz/CMMcb4IORqXt0TuvNT7k8UFBWUHktIgPWHKvZ3AZzy9CkcPnHYnyEaY4ypQUglr0f++wgiQqe4TvyU+1Pp8ZYtYfOJiv1dACt+vYLYyKB8AN0YY5qskEpezSObU1RcVKHfKyEBtlOx5nXw2EEy92ciIv4O1RhjTDVCKnnde/a9xETGVOj3ioo/zKHwLE7vcHqZ8jvzdjJ73Wx/h2mMMaYGIZW8Jrw/gW+2fVNmuDzA3ohviTtyOs3Cm5Up369dP54abVNDGWNMYxNSow3/9vO/0TqmNTvzdrLkpyWlx38qWk50TsXBGh+s/4CC4gJ7UNkYYxqZkEpeJ4pOsP/o/jJ9XkdOHGHZgY+QnRVXS+7VuhfFWuzvMI0xxtQgpJoNP1j/AV/89AU9WvVg88HNvJfxHskvJtMpviNF635RoXy/dv04o+MZAYjUGGNMdUJqMUpviU8mkhCdwPNjnmdw+xG0aQNHj5Ytc/FbFzN16FR+3vPn9X4/Y4xpbIJ5McqQSl5fb/2ajfs3csMZN7Bp/ya6texGZHgkqhAV5azvFRV1sryqoihhElIVVGNMiAjm5BVSv5UTohPoEt8FgJ6texIZHgmAiPOgcm5u2fKz1s6isLjQ32EaY4ypgavJS0RGi8gGEckUkQeqKHO1iGSIyDoRedvNePq168eFPS6s9FxCAhw8WPbY/Kz5CPaAsjHGuEHSpZWkS39Jlx6SXrsmLteaDUUkHPgRGAlsA5YDE1Q1w6tMb2A2cIGqHhCR9qq6p7r71qfZcM3uNdz+0e18cfMXFc6lpsLf/w5nVZwhyhhjmqRANBtKurQEJgMTgGbAXiAa6AB8A7yoabqwpvu4WfMaDGSqapaqngBmAeWnaP818IKqHgCoKXHVV792/fjk+k8qPVe+5rVp/yZu+/dtboZjjDGh6D1gK3CupumpmqbnaJqmapp2BR4DLpN0uaWmm7j5nFcXT4AltgFDypXpAyAiXwLhwDRVrTy7NIAwCWPxlsWM6T2mwrnyfV7tYtsx8YyJboVijDEhSdN0ZDXnVgIrfblPoB9SjgB6AyOARGCxiJyuqmV6n0RkEjAJoFmzZuXv4TNBmL5iOqN7ja4wgjAhoWzyCpMwUjqk1Pm9jDHG1EzSpR1wFxADTNc03ejLdW42G24HunrtJ3qOedsGzFXVAlXdjNNH1rv8jVR1hqqmqmpqRETd862IMHfC3EqHvrdsWbbZcPa62Tzy30fq/F7GGNMUiEiCiLwnIj+IyHoRGSYirUVkvohs9PxsVY+3+BvwKfAB4POgPTeT13Kgt4gkiUgzYDwwt1yZ/8OpdSEibXGaEbNcjIn75t/HloNbKhwv32x488CbeXzk426GYowxweAZ4BNV7QucAawHHgAWqGpvYIFn3yeSLp9KupzndagZkO15RVV2TWVcS16qWgjcgZNR1wOzVXWdiDwiIpd6in0K5IhIBrAQ+B9VzXErJoALki6gRVSLCsfLD9j494//5od9P7gZijHGNGoi0hI4D3gVQFVPeLp1LgPe8BR7A6g4v17VrgYukXR5R9KlJ/AH4M84SfK3vt7E1T4vVZ0HzCt37GGvbQV+53n5xaieoyqdbLd8zWv/0f0cKzzmr7CMMaYxSsIZyv6/InIGzmCKu4AOqrrTU2YXzjB3n2ia5gL/I+nSA3gU2AHcoWl6sPorywr0gA2/u3nuzVyYdCHXp1xf5nj5mtcNZ9zg58iMMcbvIkRkhdf+DFWd4X0eOBO4U1WXisgzlGsiVFUVEZ8fGPbUtm4HTgD3AD2BdyVdPgJe0DQt8uU+ITU9FMCrl75aIXFBxZrXqDdHsfvwbj9GZowxfldYMhjO85pR7vw2YJuqLvXsv4eTzHaLSCcAz8/aPKP7DvAvnK6if2iaLtE0HQUcBD7z9SYhl7x+2PcDy7cvr3C8fM3rryP/SuuY1n6MzBhjGhdV3QVsFZFTPYcuBDJwBt+VPAg7EZhTi9tGAZtxBmg0L32vNJ0JjPP1JiHXbPhT7k/k5OdwVpey80B517wKigpoFt6sdOJeY4wJYXcCb3lGjWcBN+FUfGaLyC3AFpxBGL76LfA8TrNhmWmMNE2PVnpFJUJqSZTq7N8PPXvCgQOw7dA2bp5zM5/9yucarDHGBJ1gXhIl5JLX97u+583Vb/LEz58oc7yw0FnLq6AAwkKuMdUYE4oCNDHvh8BLwKeapgXlzvUAbgSyNU1fq+4+Idds2LVlV67qf1WF4xER0Lw5HD4Mm46s4secH7nmtGsCEKExxjRpv8Z5POoZSZf9nJxVvjuwCXhe07TGPrSQS14J0Qn0bNWz0nMl/V6RUZHERMb4OTJjjGn6NE13AfcB90m6dAc6AUeBHzVN8329T8g1Gx48dpBBMwaxacqmCuf694d334XTTqtPhMYYExyCuc8r5Hp3EqITKk1ccHK4/H3z7+Ot1W/5OTJjjDG+CrnkBfDKt6+QdzyvwvGSZsOHzn2IS0+9tJIrjTHGNAYhmbx25O3gRNGJCsdLal4bczYiIgGIzBhjQoOkyyWSXsn6VD4KyeT18M8epk3zNhWOl9S8Xlj+AnuP7A1AZMYYEzKuATZKuvxF0qVvbS8OuQEbAPd+di/j+oxjRPcRZY4/+CDEx8NDD9Xr9iYACgoK2LZtG8eO2UoA1YmOjiYxMZHISJs9xgR+wIakSwtgAs6sHQr8L/COpmnFfp3y11aXvES4XpU3PdvDVfnS69wdqjxf3+BrqyGS1485P9I+tj0J0Qlljj/2GOw9cIxj593DC2NfqNd7GP/avHkz8fHxtGnTxpp8q6Cq5OTkkJeXR1JSUqDDMY1AoJMXgKRLG+BXwN04az/2Ap7VNH2uuutqajb0Xmer/I1urm2QjUWbmDYUFBVUOO70eSnDug4LQFSmPo4dO2aJqwYiQps2bax2ahoFSZdLJV0+ABYBkcBgTdOLcVZrvqem62t6SFmq2K5sP2i8ufpNwsPCuWPwHWWOt2wJh3OjK10yxTR+lrhqZv9GphG5AnhK03Sx90FN03xJl1tqurimmpdWsV3ZftC4a+hdFRIXOMlrIx9xxewrAhCVCXYiwvXXn/zDp7CwkHbt2jFunLPKw+7duxk3bhxnnHEGycnJjBkzBoDs7GxiYmIYMGBA6WvmzJkB+QzG+NE0YFnJjqRLjGfGDTRNF9R0cU01r74irMapZfX0bOPZ71GHYBuFjL0ZfL/reyacPqHM8YQEaLZlLG/98qIARWaCWWxsLGvXruXo0aPExMQwf/58unTpUnr+4YcfZuTIkdx1110ArF69uvRcz549+e677/weszEB9E/gbK/9Is+xsyovXlZNNa9+wCU4C4SVbJfsJ9c20sYiTMIqbT5p2RL2sIbsg9n+D8o0CWPGjOGjjz4C4J133mHChJN/IO3cuZPExMTS/ZSUFL/HZ0wjEqFpWvrArWe7ma8XV5u8VNni/QIO4ywB3dazXy0RGS0iG0QkU0QeqOT8jSKyV0S+87xu9TXw+ujbti/jTxtf5tgP+35gVd4nHIhcw/q96/0RhmmCxo8fz6xZszh27BirV69myJAhpecmT57MLbfcwvnnn8+jjz7Kjh07Ss9t2rSpTLPhkiVLAhG+Mf60V9KldCojSZfLgH2+Xlxts6EI/wYeUGWtCJ2Ab4EVOE2IM1R5uuprJRx4ARgJbAOWi8hcVc0oV/RdVa3YAeWirblbue5f17H4psWs2b2GguICioqLOBq2lxMrf8Xl/fwZjXGDG+MSfHkkMiUlhezsbN55553SPq0So0aNIisri08++YSPP/6YgQMHsnbtWsCaDU1Iug14S9LleZyuqK3ADb5eXFOfV5Iqaz3bNwHzVblBhHjgS6g6eQGDgUxVzQIQkVnAZUD55OV3HeM68o/L/wHAT7k/cfjEYa457RpSO5/F7cedBSntGc7gFshn7y+99FLuvfdeFi1aRE5OTplzrVu35tprr+Xaa69l3LhxLF68mEGDBgUoUmMCR9N0EzBU0iXOs3+4NtfXlLy8H4a6EHgZQJU8EYpruLYLTiYtsQ0YUkm5K0TkPOBHYKqqbi1fQEQmAZMAmjXzuUm0SpHhkew6vIu9+XsZ22es1/ucnCKqbdt6v40JUTfffDMJCQmcfvrpLFq0qPT4559/ztChQ2nevDl5eXls2rSJbt26BS5QYwJM0mUs0B+IlnSnuUTT9BFfrq1pwMZWEe4U4XKcvq5PAESIwXmorL4+BLqragowH3ijskKqOkNVU1U1NSKiYdbPvPvTu7lj3h3kF5Rd+6xkcl5j6ioxMZEpU6ZUOL5y5UpSU1NJSUlh2LBh3HrrrZx1ljOwqnyf17PPPuvvsI3xK0mX6TjzG96J02x4FXCKz9fXMD1Ue+ARnJUuX1DlM8/x84FBqvy16mtlGDBNVUd59h8EUNU/V1E+HNivqi2rC7ghpoeqzplnwssvg7XkBJf169fTr591VvrC/q1MiUBODyXpslrTNMXrZxzwsabpub5cX201RpU9OJ1q5Y8vBBbWcO/lQG8RSQK2A+OBa8sEL9JJVXd6di/FmdcqoKzmZYwxflEyT1m+pEtnIAenouSTmkYbzq3uvCpVrtioqoUicgfwKRAOvKaq60TkEWCFqs4FpojIpUAhsB+40dfA3VLS52WMMcZVH0q6JABP4IxkVzzjKnxRUwfSMJxBF+8AS6nlfIaqOg+YV+7Yw17bDwIP1uaebmvVCvb5/KSBMcaY2vIsQrlA0/Qg8L6ky7+BaE1Tn6sONQ3Y6Ag8BJwGPIPzzNY+Vf6ryn/rGHejduaZsGxZzeWMMcbUjaZpMc5zwCX7x2uTuKDmGTaKVPlElYnAUCATWCSCXx8q9qfzzoPFi2suZ4wxpl4WSLpcIel1m1KgxojEjm8AABlQSURBVHHnIkQBY3FWu+wOPAt8UJc3CwanneY0G+7aBR07BjoaY4xpsn6Ds2ZkoaTLMZxuKdU0beHLxTUN2JiJ02Q4D0j3mm2jyQoLg+HDYckSuOqqQEdjgkVOTg4XXnghALt27SI8PJx27doBsGzZsmofrl+xYgUzZ86s8dmus88+m6+++qrhgjYmgDRN4+tzfU3PeRUDJQ9VeRd0MqTiU4ZsSG4/5wXwxBPw00/wXLWLUJvGpDE9uzRt2jTi4uK49957S48VFhbSUA/Y11dj+rcygRXg57zOq+x4+cUpq1JTn1eYKvGeVwuvV3wgEpe/WL+XaQg33ngjt912G0OGDOG+++5j2bJlDBs2jIEDB3L22WezYcMGABYtWlS6YOW0adO4+eabGTFiBD169ChTG4uLiystP2LECK688kr69u3LddddR8kfofPmzaNv374MGjSIKVOmlN7XmEbof7xef8CZcWmarxc3jj8FG5kzz4SsLDhwwBk6b0xdbdu2ja+++orw8HAOHTrEkiVLiIiI4D//+Q8PPfQQ77//foVrfvjhBxYuXEheXh6nnnoqt99+O5HlZopetWoV69ato3PnzgwfPpwvv/yS1NRUfvOb37B48WKSkpLKrCVmTGOjaXqJ976kS1eqn+y9DEtelYiMhKFD4csvwf5wDU4lk3w2JE2r/VT1V111FeHh4QDk5uYyceJENm7ciIhQUFBQ6TVjx44lKiqKqKgo2rdvz+7du8ssYgkwePDg0mMDBgwgOzubuLg4evToQVJSEgATJkxgxowZtY7ZmADZhrPosU8seVWhpOnQkldwqkuicUNs7MnuhD/84Q+cf/75fPDBB2RnZzNixIhKr4mKiirdDg8Pp7CwsE5ljGnMJF2e4+RYijBgAM5MGz6x5FWFc8+FByqs/WxM3eXm5tKlSxcAXn/99Qa//6mnnkpWVhbZ2dl0796dd999t8Hfw5gGtMJruxB4R9P0S18vrmmGjZA1ZAisWQMuD2w0IeS+++7jwQcfZODAga7UlGJiYnjxxRcZPXo0gwYNIj4+npYtq12kwZhAeg94U9P0DU3Tt4BvJF2a+3pxtUPlGyN/DJUvcc45kJ4Onsd3TCNmw78dhw8fJi4uDlVl8uTJ9O7dm6lTp5YpY/9WpkSAh8p/A1xUsoKyZ0mUzzRNz/bleqt5VeOcc5xBG8YEi5dffpkBAwbQv39/cnNz+c1vfhPokIypSnRJ4gLwbPtc87LkVY1+/SAzM9BRGOO7qVOn8t1335GRkcFbb71F8+Y+/y4wxt+OSLqcWbIj6TIIOOrrxTZgoxo9ejirKhtjTKjyrHK/AtiuquM8CwzPAtoAK4FfqeqJOtz6buCfki47cGZt6ghc4+vFVvOqRo8ezsPKxhgTwu6i7Cr3jwNPqWov4ABwS11uqmm6HOgL3A7cBvTTNF3p6/WWvKrRqZMzy0Z+fqAjMcYY/xORRJxVRV7x7AtwAc5IQYA3gF/U6d7pMhmI1TRdq2m6FoiTdPmtr9db8qpGWBh07w6bNwc6EmOMCYingfuAYs9+G+CgqpY867EN6FLHe//as5IyAJqmB4Bf+3qxJa8aWNOh8UVOTg4DBgxgwIABdOzYkS5dupTunzhRc3fAokWLyix3Mn36dGbOnOlmyMYARIjICq/XpJITIjIO2KPqe1NeLYV7L0Qp6RIOVL12UDmuDtgQkdHAM0A48IqqPlZFuStwqqFnqeqKysoEiiUv44s2bdrw3XffAZUviVKTRYsWERcXx9lnO4+43Hbbba7EaUw5haqaWsW54cClIjIGiAZa4Pw+TxCRCE/tKxHYXsf3/gR4V9LlJc/+bzzHfOJazcszQuUF4GIgGZggIsmVlIvH6RBc6lYs9WHJy9TVypUr+dnPfsagQYMYNWoUO3fuBODZZ58lOTmZlJQUxo8fT3Z2NtOnT+epp55iwIABLFmyhGnTpvHXv/4VgBEjRnD//fczePBg+vTpw5IlSwDIz8/n6quvJjk5mcsvv5whQ4awYkWj+tvPBDFVfVBVE1W1OzAe+FxVrwMWAld6ik0E5tTxLe4HPscZsHE7sABneRSfuFnzGgxkqmoWgIjMAi4DMsqV+yPO6BWfg/anHj3g888DHYUJNqrKnXfeyZw5c2jXrh3vvvsuv//973nttdd47LHH2Lx5M1FRURw8eJCEhARuu+22MrW1BQsWlLlfYWEhy5YtY968eaSnp/Of//yHF198kVatWpGRkcHatWsZMGBAID6qCT33A7NE5E/AKuDVutxE07QYmO55IelyLvAcMNmX693s8+oCbPXar9CxJyJnAl1V9SMX46gXq3kFp2mLpjFt0TQA+jzXhx9zfmTljpUMmjEIgHs+vYe/ffU3ADr/rTM78nawKHsRI14fAcCkDycxY6WznEj8n+PJO55Xq/c/fvw4a9euZeTIkQwYMIA//elPbNu2DYCUlBSuu+463nzzTZ9XV/7lL38JwKBBg8jOzgbgiy++YPz48QCcdtpppKSk1CpGY3ylqotUdZxnO0tVB6tqL1W9SlWP1/W+ki4DJV3+IumSDTwC/ODrtQF7SFlEwoAngRt9KDsJmATQrJnP/XkNIinJGW2oCtLwS0QZl0wbMa10+8c7fyzdXjnJ6Xv+26i/lR7bcc8OADrHd2bRjYsAmHHJyXWw8h6sXeICp+bVv39/vv766wrnPvroIxYvXsyHH37Io48+ypo1a2q8X8kSKLb8iQl2ki59gAme1z7gXUA0Tc+vzX3crHltB7p67Zfv2IsHTgMWiUg2MBSYKyIVOg9VdYaqpqpqqq9/qTaUuDiIj4ddu/z6tibIRUVFsXfv3tLkVVBQwLp16yguLmbr1q2cf/75PP744+Tm5nL48GHi4+PJy6tdkhw+fDizZ88GICMjw6ckaEwj8APOs2LjNE3P0TR9Diiq7U3czATLgd6eqUS243T4XVtyUlVzgbYl+yKyCLi3sY02hJNNh506BToSEyzCwsJ47733mDJlCrm5uRQWFnL33XfTp08frr/+enJzc1FVpkyZQkJCApdccglXXnklc+bM4bnnnvPpPX77298yceJEkpOT6du3L/3797clUEww+CVOPlgo6fIJzlRTtW7XcnVJFM8Qy6dxhsq/pqqPisgjwApVnVuu7CJ8SF7+XBKlxPXXw6hR8Ktf+fVtTS2E4jIfRUVFFBQUEB0dzaZNm7jooovYsGFDjU3rofhvZSoX4CVRYnEG8U3AqYnNBD7QNP3Ml+tdbYNT1XnAvHLHHq6i7Ag3Y6kPG7RhGqP8/HzOP/98CgoKUFVefPFFv/cJG1NXmqZHgLeBtyVdWgFX4YxkDHzyaip69ICFCwMdhTFlxcfH23NdpknwTA01w/PyiU0P5QOreRljTONiycsHlryCg5v9t02F/RuZpsKSlw86d4b9++Goz2t8Gn+Ljo4mJyfHfjlXQ1XJyckhOjo60KEYU2/W5+WDsDA45RTIzgYbpNU4JSYmsm3bNvbu3RvoUBq16OhoEhMTAx2GMfVmyctHJU2Hlrwap8jISJKSkgIdhjHGT6zZ0Ec9esCmTYGOwhhjDFjy8pklL2OMaTwsefmoZ08bcWiMMY2FJS8f9eoFmZmBjsIYYwxY8vJZjx7OaMOiWs99bIwxpqFZ8vJRTAy0bQue9QSNMcYEkCWvWrCmQ2OMaRwsedWCJS9jjGkcLHnVgiUvY4xpHCx51YIlL2OMaRwsedWCJS9jjGkcLHnVQs+eziwbxcWBjsQYY0KbJa9aiIuDli1h585AR2KMMaHNklctWdOhMcYEnqvJS0RGi8gGEckUkQcqOX+biKwRke9E5AsRSXYznoZgycsYYwLPteQlIuHAC8DFQDIwoZLk9Laqnq6qA4C/AE+6FU9DseRljDGB52bNazCQqapZqnoCmAVc5l1AVQ957cYCjX4Nd0texhgTeG6upNwF2Oq1vw0YUr6QiEwGfgc0Ay6o7EYiMgmYBNCsWbMGD7Q2LHkZY0zgBXzAhqq+oKo9gfuB/6+KMjNUNVVVUyMi3My3NevZ00le2ujriMYY03S5mby2A1299hM9x6oyC/iFi/E0iIQEiI6GPXsCHYkxxoQuN5PXcqC3iCSJSDNgPDDXu4CI9PbaHQtsdDGeBmNNh8YYE1iuJS9VLQTuAD4F1gOzVXWdiDwiIpd6it0hIutE5Ducfq+JbsXTkCx5GWNMYLnagaSq84B55Y497LV9l5vv7xZLXsYYE1gBH7ARjCx5GWNMYFnyqoO+fWHNmkBHYYwxocuSVx2kpEB2NuTmBjoSY4wJTZa86iAyEgYNgqVLAx2JMcaEJktedXT22fDVV4GOwhhj3CMiXUVkoYhkeEaG3+U53lpE5ovIRs/PVv6OzZJXHVnyMsaEgELgHlVNBoYCkz0TrD8ALFDV3sACz75fWfKqo6FDnWbDoqJAR2KMMe5Q1Z2q+q1nOw/nmd0uOJOsv+Ep9gYBmB3JklcdtWsHHTpARkagIzHGmDqLEJEVXq9JVRUUke7AQGAp0EFVS9aU3wV0cD3ScgI7y22QO/ts+PprOP30QEdijDF1UqiqqTUVEpE44H3gblU9JCKl51RVRcTvU5Vbzasehg2zfi9jTNMmIpE4iestVf2X5/BuEenkOd8J8PtU5Za86sEGbRhjmjJxqlivAutV1Xul+7mcnIt2IjDH77FpkC1MFRsbq0eOHAl0GIAzWKN1a9i0Cdq2DXQ0xhhTOyKSr6qx1Zw/B1gCrAGKPYcfwun3mg10A7YAV6vqfpfDLcP6vOohPByGDHH6vS65JNDRGGNMw1LVLwCp4vSF/oylPGs2rCdrOjTGGP+z5FVPw4bBl18GOgpjjAktlrzqafhwZ5LeDz8MdCTGGBM6LHnVU1wcvPMO3HorbN0a6GiMMSY0WPJqAMOHw9SpMH48FBQEOhpjjGn6LHk1kPvug/h4SEsLdCTGGNP0uZq8RGS0iGwQkUwRqTDrsIj8zjPV/moRWSAip7gZj5vCwmDmTHj9dWfovDHGGPe4lrxEJBx4AbgYSAYmeKbS97YKSFXVFOA94C9uxeMP7dvDk0/CpEnWfGiMMW5ys+Y1GMhU1SxVPQHMwplGv5SqLlTVfM/uN0Cii/H4xTXXQGKik8SMMca4w83k1QXwHn+3zXOsKrcAH7sYj1+IwIsvwhNPQFZWoKMxxpimqVEM2BCR64FU4Ikqzk8qWW+msLDQv8HVQVKSM4Dj9tshyKaONMaYoOBm8toOdPXaT/QcK0NELgJ+D1yqqscru5GqzlDVVFVNjYgIjukYp06FQ4ecBGarLRtjTMNyM3ktB3qLSJKINAPG40yjX0pEBgIv4SQuv68H46bISPjsM9i40ekHO15pWjbGGFMXriUvVS0E7gA+BdYDs1V1nYg8IiKXeoo9AcQB/xSR70RkbhW3C0rx8TBvntMPNmYM7N0b6IiMMaZpsPW8/KCoCB56CF59Fe6+22lSjK1yBR1jjPGPmtbzaswaxYCNpi48HB5/HJYuhXXroHdvWLIk0FEZY0zwsppXAHz0kfMg87ffQocOgY7GGBOqrOZlamXsWLjpJvjVr6C4uObyxhhjyrLkFSDTpsGxY/DnPwc6EmOMCT7B8dBUExQRAW+/DampMGQIXHRRoCMyxpjgYTWvAEpMhFmz4NprnRnpjTHG+MZqXgE2YgQsWgTjxsGGDfDHPzrLqxhjjKma/ZpsBJKTnWH0//2vk8z++U84cSLQURljTONlQ+UbkRMn4IMPYPp0+OEHpznxoovgnHOc2TqMMaYh2VB50yCaNXPmQVy4ED7/HFq0cB5u7tQJRo50mhWNMcZY8mq0+vWDtDSnP2zvXrjsMhg+HJ591p4NM8YYazYMIhs3wsSJTvIaOBDatnVm6Bg5Ek49NdDRGWOCTTA3G1ryCjJFRfDhh7B9O+TkwNatznRTbdvClVfCoEHOYpjdu0Pz5oGO1hjTmFny8qNQT16VKS6Gr75yBntkZEBWFmzZAnFx0KULdO7s1NDat3debdpAy5bOq0ULJ8mVf4WHB/pTGWPcZsnLjyx5+aa4GPbtgx07nFra7t1O39mePU6NLTcXDh50Vns+ehTy8+HIEWf76FFnBpAWLU4muOhoiIpyfkZGnnxFRDiJruRnWNjJV/n9sDBnbbPqXnBy27t8Ce8y3j8r2/alTE3b5flSrqGO16Wcr/dy6/rG8h6NkVuf+5Zb6v7HpiUvP7Lk5T5VZ97FvDwnyR065OwfP+78LCg4+SoqgsJC51VcfPJVVFR2W/Xkvmrlr5L39n55D07xLuP9s7JtX8rUtF3Zv0tN5RrqeF3K1fd/ZX/8KgiyXzcNxs3P/eKLzh+PdRHMyctm2DAViEBMjPNq3z7Q0RhjTEU2VN4YY0zQseRljDEm6FjyMsYYE3RcTV4iMlpENohIpog8UMn580TkWxEpFJEr3YzFGGNM0+Fa8hKRcOAF4GIgGZggIsnliv0E3Ai87VYcxhhjmh43a16DgUxVzVLVE8As4DLvAqqaraqrAZutzxhjGqGaWtACxc3k1QXY6rW/zXOs1kRkkoisEJEVhYWFDRKcMcaY6vnYghYQQTFgQ1VnqGqqqqZG1PVpPGOMMbVVYwtaoLiZvLYDXb32Ez3HjDHGBIcGa0FraG5WY5YDvUUkCSdpjQeure9N8/PzVUSO1uMWEUAotj3a5w49ofrZQ/VzQ+0/e4yIrPDan6GqMxo4Jle4lrxUtVBE7gA+BcKB11R1nYg8AqxQ1bkichbwAdAKuERE0lW1fw33rVdtUURWqGpqfe4RjOxzh55Q/eyh+rnBlc/eaFvQXO1AUtV5wLxyxx722l6O849hjDGm8XGlBa0h2OgHY4wxlaqqBS3AYQGhmbyCoj3XBfa5Q0+ofvZQ/dzgwmevrAWtMQi69byMMcaYoHjOyxhjjPEWMsmrsU5x0tBEpKuILBSRDBFZJyJ3eY63FpH5IrLR87NVoGN1i4iEi8gqEfm3Zz9JRJZ6vvt3RaRZoGNsaCKSICLvicgPIrJeRIaFyncuIlM9/62vFZF3RCS6KX7nIvKaiOwRkbVexyr9jsXxrOfzrxaRMwMXuTtCInk15ilOXFAI3KOqycBQYLLnsz4ALFDV3sACz35TdRew3mv/ceApVe0FHABuCUhU7noG+ERV+wJn4Hz+Jv+di0gXYAqQqqqn4QwqGE/T/M5fB0aXO1bVd3wx0NvzmgT83U8x+k1IJC8a8RQnDU1Vd6rqt57tPJxfYl1wPu8bnmJvAL8ITITuEpFEYCzwimdfgAuA9zxFmtxnF5GWwHnAqwCqekJVDxIi3znOwLMYEYkAmgM7aYLfuaouBvaXO1zVd3wZMFMd3wAJItLJP5H6R6gkr0Y7xYmbRKQ7MBBYCnRQ1Z2eU7uADgEKy21PA/dxcqWCNsBBVS2ZdaApfvdJwF7gfz3Npa+ISCwh8J2r6nbgrzjLK+0EcoGVNP3vvERV33GT/50XKskr5IhIHPA+cLeqHvI+p84Q0yY3zFRExgF7VHVloGPxswjgTODvqjoQOEK5JsIm/J23wqllJAGdgVgqNq2FhKb6HVclVJJXo53ixA0iEomTuN5S1X95Du8uaTbw/NwTqPhcNBy4VESycZqGL8DpC0rwNClB0/zutwHbVHWpZ/89nGQWCt/5RcBmVd2rqgXAv3D+O2jq33mJqr7jJv87L1SSV+kUJ55RR+OBuQGOyRWePp5XgfWq+qTXqbnARM/2RGCOv2Nzm6o+qKqJqtod5zv+XFWvAxYCV3qKNbnPrqq7gK0icqrn0IVABiHwneM0Fw4Vkeae//ZLPnuT/s69VPUdzwVu8Iw6HArkejUvNgkh85CyiIzB6Q8pmeLk0QCH5AoROQdYAqzhZL/PQzj9XrOBbsAW4GpVLd/522SIyAjgXlUdJyI9cGpirYFVwPWqejyQ8TU0ERmAM0ilGZAF3ITzx2mT/85FJB24Bmek7SrgVpz+nSb1nYvIO8AIoC2wG0gD/o9KvmNPIn8epwk1H7hJVVdUdt9gFTLJyxhjTNMRKs2GxhhjmhBLXsYYY4KOJS9jjDFBx5KXMcaYoGPJyxhjTNCx5GVMNUSkSES+83o12OS2ItLde4ZwY4zvQnElZWNq46iqDgh0EMaYsqzmZUwdiEi2iPxFRNaIyDIR6eU53l1EPvesobRARLp5jncQkQ9E5HvP62zPrcJF5GXPelSfiUhMwD6UMUHEkpcx1Ysp12x4jde5XFU9HWcmg6c9x54D3lDVFOAt4FnP8WeB/6rqGTjzDq7zHO8NvKCq/YGDwBUufx5jmgSbYcOYaojIYVWNq+R4NnCBqmZ5JkLepaptRGQf0ElVCzzHd6pqWxHZCyR6T1HkWbJmvmchQUTkfiBSVf/k/iczJrhZzcuYutMqtmvDe769Iqwf2hifWPIypu6u8fr5tWf7K5wZ7QGuw5kkGZwl2m8HEJFwz+rHxpg6sr/yjKlejIh857X/iaqWDJdvJSKrcWpPEzzH7sRZ0fh/cFY3vslz/C5ghojcglPDuh1n5V9jTB1Yn5cxdeDp80pV1X2BjsWYUGTNhsYYY4KO1byMMcYEHat5GWOMCTqWvIwxxgQdS17GGGOCjiUvY4wxQceSlzHGmKBjycsYY0zQ+f8Be5Sioc52D34AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_data():\n",
    "    x = range(0, EPOCHS)\n",
    "    fig, ax2 = plt.subplots()\n",
    "    ax2.set_xlabel('Epoch')\n",
    "    ax2.set_ylabel('MSE', color='blue')\n",
    "    line, = ax2.plot(x, MSE, '-', c='blue', lw='1', label='MSE')\n",
    "    ax1 = ax2.twinx()\n",
    "    ax1.set_ylabel('Accuracy (%)', color='green')\n",
    "    line2, = ax1.plot(x, TRP, '-', c='green', lw='1', label='Training')\n",
    "    line3, = ax1.plot(x, TEP, ':', c='green', lw='1', label='Testing')\n",
    "    fig.tight_layout()\n",
    "    fig.legend(loc='center')\n",
    "    ax1.set_ylim(0, 101)\n",
    "    plt.show()\n",
    "    plt.clf()\n",
    "    \n",
    "plot_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see over time the mean squared error is diminished, just as it would were you using backpropagation; as the mean squared error is inversely proportional to the accuracy of the network, you also see accuracy rise.\n",
    "\n",
    "## Final Thoughts\n",
    "\n",
    "It's interesting to think about a metaheuristic like particle swarm optimization algorithms in training a network. A lot of literature and even what you may have learned in school were you to take a computer science degree use PSO in different ways. At least for me, using one to train a neural network isn't standard curriculum but it still has merit to its study.\n",
    "\n",
    "Using this method to train a neural network, we can achieve mostly comparable results. What I didn't mention in this tutorial is how slow it can be compared to backpropagation. However, depending on parameter choices, it tends to converge faster but not necessarily to the same found optimum. It would be useful to define some termination condition (for example, terminate if MSE is below a certain threshold) as often you do not need to run for many epochs but this is problem specific.\n",
    "\n",
    "Further consideration should be given to the activation function used. ReLU has many pitfalls in backpropagation trained networks but performs great in metaheuristic trained networks. I chose it specifically because this is one scenario where there are few drawbacks to doing so. You should still look into other types of activation functions as well as experiment with differing PSO parameters to potential get better results.\n",
    "\n",
    "The difficulty in using a metaheuristic like PSO is that the search space isn't known a priori. We can make some educated assumptions that idealized network weights are relatively small in magnitude and base some of our conditions around that: for example, initializing small and not updating fitness of a particle if it goes too far from the origin. PSO cares a little less about search boundaries than genetic algorithms but it's still an important thing to consider. This is why I've also normalized the data: doing this means weights remain closer to the origin, thus less orphan cases to catch were network weights large in comparison.\n",
    "\n",
    "## Finished Code\n",
    "\n",
    "Here is the full uncommented code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\n",
    "import pandas as pd\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "class Particle:\n",
    "    def __init__(self, pos, vel):\n",
    "        self.pos, self.vel = pos, vel\n",
    "        network = initialize_network(self.pos)\n",
    "        self.fit = mse(network)\n",
    "        self.best_pos, self.best_fit = self.pos, self.fit\n",
    "\n",
    "    def set_pos(self, pos):\n",
    "        self.pos = pos\n",
    "        if not any(p < -7.00 for p in pos)\\\n",
    "        and not any(p > 7.00 for p in pos):\n",
    "            network = initialize_network(self.pos)\n",
    "            fitness = mse(network)\n",
    "            if fitness < self.best_fit:\n",
    "                self.fit = fitness\n",
    "                self.best_fit = self.fit\n",
    "                self.best_pos = self.pos\n",
    "\n",
    "    def set_vel(self, vel):\n",
    "        self.vel = vel\n",
    "\n",
    "    def get_pos(self):\n",
    "        return self.pos\n",
    "\n",
    "    def get_vel(self):\n",
    "        return self.vel\n",
    "\n",
    "    def get_best_pos(self):\n",
    "        return self.best_pos\n",
    "\n",
    "    def get_fit(self):\n",
    "        return self.fit\n",
    "\n",
    "def pso(dim, epochs, swarm_size, ic, cc, sc):\n",
    "    swarm = initialize_swarm(swarm_size, dim)\n",
    "    for e in range(1, epochs+1):\n",
    "        swarm_best = get_swarm_best(swarm)\n",
    "        MSE.append(swarm_best[0])\n",
    "        TRP.append(performance_measure(swarm_best[1], TRAIN))\n",
    "        TEP.append(performance_measure(swarm_best[1], TEST))\n",
    "        move_particles(swarm, dim, ic, cc, sc)\n",
    "    return get_swarm_best(swarm)[0]\n",
    "\n",
    "def move_particles(swarm, dim, ic, cc, sc):\n",
    "    swarm_best = get_swarm_best(swarm)\n",
    "    for particle in swarm:\n",
    "        new_pos = [0 for _ in range(dim)]\n",
    "        new_vel = [0 for _ in range(dim)]\n",
    "        for d in range(dim):\n",
    "            r_1 = random.uniform(0.00, 1.00)\n",
    "            r_2 = random.uniform(0.00, 1.00)\n",
    "            weight = ic * particle.get_vel()[d]\n",
    "            cognitive = cc * r_1\n",
    "            cognitive *= (particle.get_best_pos()[d] - particle.get_pos()[d])\n",
    "            social = sc * r_2\n",
    "            social *= (swarm_best[1][d] - particle.get_pos()[d])\n",
    "            new_vel[d] = weight + cognitive + social\n",
    "            new_pos[d] = particle.get_pos()[d] + new_vel[d]\n",
    "        particle.set_pos(new_pos)\n",
    "        particle.set_vel(new_vel)\n",
    "\n",
    "def feed_forward(network, example):\n",
    "    layer_input, layer_output = example, []\n",
    "    for layer in network:\n",
    "        for neuron in layer:\n",
    "            summ = summing_function(neuron, layer_input)\n",
    "            layer_output.append(activation_function(summ))\n",
    "        layer_input, layer_output = layer_output, []\n",
    "    return layer_input\n",
    "\n",
    "def summing_function(weights, inputs):\n",
    "    bias = weights[-1]\n",
    "    summ = 0.00\n",
    "    for i in range(len(weights)-1):\n",
    "        summ += (weights[i] * float(inputs[i]))\n",
    "    return summ + bias\n",
    "\n",
    "def activation_function(z):\n",
    "    return z if z >= 0 else 0\n",
    "\n",
    "def sse(actual, target):\n",
    "    summ = 0.00\n",
    "    for i in range(len(actual)):\n",
    "        summ += (actual[i] - target[i])**2\n",
    "    return summ\n",
    "\n",
    "def mse(network):\n",
    "    training = TRAIN\n",
    "    summ = 0.00\n",
    "    for example in training:\n",
    "        target = [0 for _ in range(CLASSES)]\n",
    "        target[int(example[-1])] = 1\n",
    "        actual = feed_forward(network, example)\n",
    "        summ += sse(actual, target)\n",
    "    return summ / len(training)\n",
    "\n",
    "def performance_measure(particle, data):\n",
    "    network = initialize_network(particle)\n",
    "    correct, total = 0, 0\n",
    "    for example in data:\n",
    "        if check_output(network, example) == float(example[-1]):\n",
    "            correct += 1\n",
    "        total += 1\n",
    "    return 100*(correct / total)\n",
    "\n",
    "def check_output(network, example):\n",
    "    output = feed_forward(network, example)\n",
    "    return output.index(max(output))\n",
    "\n",
    "def initialize_network(p):\n",
    "    n, h, o = FEATURES, HIDDEN_SIZE, CLASSES\n",
    "    part = iter(p)\n",
    "    neural_network = []\n",
    "    neural_network.append([[next(part) for i in range(n+1)] for j in range(h)])\n",
    "    neural_network.append([[next(part) for i in range(h+1)] for j in range(o)])\n",
    "    return neural_network\n",
    "\n",
    "def initialize_swarm(size, dim):\n",
    "    swarm = []\n",
    "    for _ in range(size):\n",
    "        position = [random.uniform(-1.00, 1.00) for _ in range(dim)]\n",
    "        velocity = [0 for _ in range(dim)]\n",
    "        particle = Particle(position, velocity)\n",
    "        swarm.append(particle)\n",
    "    return swarm\n",
    "\n",
    "def get_swarm_best(swarm):\n",
    "    best_fit = swarm[0].get_fit()\n",
    "    best_pos = swarm[0].get_pos()\n",
    "    for particle in swarm:\n",
    "        if particle.get_fit() < best_fit:\n",
    "            best_fit = particle.get_fit()\n",
    "            best_pos = particle.get_pos()\n",
    "    return best_fit, best_pos\n",
    "\n",
    "def load_data(filename):\n",
    "    df = pd.read_csv(filename, header=None, dtype=float)\n",
    "    for features in range(len(df.columns)-1):\n",
    "        df[features] = (df[features] - df[features].mean())/df[features].std()\n",
    "    train = df.sample(frac=0.70).fillna(0.00)\n",
    "    test = df.drop(train.index).fillna(0.00)\n",
    "    return train.values.tolist(), test.values.tolist()\n",
    "\n",
    "def plot_data():\n",
    "    x = range(0, EPOCHS)\n",
    "    fig, ax2 = plt.subplots()\n",
    "    ax2.set_xlabel('Epoch')\n",
    "    ax2.set_ylabel('MSE', color='blue')\n",
    "    line, = ax2.plot(x, MSE, '-', c='blue', lw='1', label='MSE')\n",
    "    ax1 = ax2.twinx()\n",
    "    ax1.set_ylabel('Accuracy (%)', color='green')\n",
    "    line2, = ax1.plot(x, TRP, '-', c='green', lw='1', label='Training')\n",
    "    line3, = ax1.plot(x, TEP, ':', c='green', lw='1', label='Testing')\n",
    "    fig.tight_layout()\n",
    "    fig.legend(loc='center')\n",
    "    ax1.set_ylim(0, 101)\n",
    "    plt.show()\n",
    "    plt.clf()\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    MSE = []\n",
    "    TRP = []\n",
    "    TEP = []\n",
    "    W = 0.3\n",
    "    C_1 = 1.5\n",
    "    C_2 = 1.2\n",
    "    TRAIN, TEST = load_data('../data/wine.csv')\n",
    "    FEATURES = len(TRAIN[0][:-1])\n",
    "    CLASSES = len(list(set([c[-1] for c in (TRAIN+TEST)])))\n",
    "    HIDDEN_SIZE = 5\n",
    "    DIMENSIONS = (HIDDEN_SIZE * (FEATURES+1)) + \\\n",
    "        (CLASSES * (HIDDEN_SIZE+1))\n",
    "    SWARM_SIZE = 100\n",
    "    EPOCHS = 100\n",
    "    NETWORK = pso(DIMENSIONS, EPOCHS, SWARM_SIZE, W, C_1, C_2)\n",
    "    plot_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may wish to catch for errors (as I have done none here), and perhaps accept some of the parameters as command line arguments. This code, as-is, should run in a terminal but allows for zero customization unless you change the code manually. Some parameters might include hidden layer size, PSO coefficients, how many epochs, etc. This is just a basic implementation to springboard off of and experiment with.\n",
    "\n",
    "You can save this code as `particle_network.py` and execute it as one of the below:\n",
    "\n",
    "```\n",
    "$ python3 particle_network.py\n",
    "$ ./particle_network.py\n",
    "```\n",
    "\n",
    "You can also find the code along with a commented version in `./code`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
